Merge pull request #35 from pjpegion/ocn_stoch

Ocean stochastic physics and netCDF restarts

CMakeLists.txt

@@ -16,11 +16,8 @@ list(APPEND _stoch_phys_srcs
   mpi_wrapper.F90
   halo_exchange.fv3.F90
   plumes.f90
-
-  stochy_gg_def.f
-  stochy_resol_def.f
-  stochy_layout_lag.f
   four_to_grid_stochy.F
+  fftpack_stochy.f
   glats_stochy.f
   sumfln_stochy.f
   gozrineo_stochy.f
@@ -35,7 +32,6 @@ list(APPEND _stoch_phys_srcs
   dozeuv_stochy.f
   dezouv_stochy.f
   mersenne_twister.F
-
   spectral_layout.F90
   getcon_spectral.F90
   stochy_namelist_def.F90
@@ -53,20 +49,25 @@ list(APPEND _stoch_phys_srcs
 )
 
 add_library(stochastic_physics ${_stoch_phys_srcs})
-if(32BIT)
-   add_dependencies(stochastic_physics FMS::fms_r4)
-else()
-   add_dependencies(stochastic_physics FMS::fms_r8)
-endif()
+add_dependencies(stochastic_physics fms)
 
 target_compile_definitions(stochastic_physics PRIVATE INTERNAL_FILE_NML)
 
-target_link_libraries(stochastic_physics PUBLIC sp::sp_d)
-if(32BIT)
-   target_link_libraries(stochastic_physics PUBLIC FMS::fms_r4)
-else()
-   target_link_libraries(stochastic_physics PUBLIC FMS::fms_r8)
-endif()
+target_link_libraries(stochastic_physics PUBLIC fms)
 if(OpenMP_Fortran_FOUND)
-  target_link_libraries(stochastic_physics PUBLIC OpenMP::OpenMP_Fortran)
+  target_link_libraries(stochastic_physics PRIVATE OpenMP::OpenMP_Fortran)
 endif()
+
+###############################################################################
+### Install
+###############################################################################
+install(
+  TARGETS stochastic_physics
+  EXPORT  stochastic_physics-config
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib)
+
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/mod DESTINATION ${CMAKE_INSTALL_PREFIX})
+
+install(EXPORT stochastic_physics-config
+    DESTINATION lib/cmake)

compns_stochy.F90

@@ -57,10 +57,12 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
 !
       namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, &
       iseed_shum,iseed_sppt,shum,shum_tau,&
-      shum_lscale,fhstoch,stochini,skeb_varspect_opt,sppt_sfclimit, &
+      shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, &
       skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, &
       skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,&
-      shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale
+      shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, &
+      epbl,epbl_lscale,epbl_tau,iseed_epbl,                                    &
+      ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt
       namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & 
       lndp_tau,lndp_lscale 
 
@@ -130,27 +132,25 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
 ! length scale.
       skeb_varspect_opt = 0
       sppt_logit        = .false. ! logit transform for sppt to bounded interval [-1,+1]
-      fhstoch           = -999.0  ! forecast interval (in hours) to dump random patterns
       stochini          = .false. ! true= read in pattern, false=initialize from seed
 
 #ifdef INTERNAL_FILE_NML
       read(input_nml_file, nml=nam_stochy)
 #else
       rewind (nlunit)
-      open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios)
+      open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios)
       read(nlunit,nam_stochy)
 #endif
 #ifdef INTERNAL_FILE_NML
       read(input_nml_file, nml=nam_sfcperts)
 #else
       rewind (nlunit)
-      open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios)
+      open (unit=nlunit, file=fn_nml, action='READ', status='OLD', iostat=ios)
       read(nlunit,nam_sfcperts)
 #endif
 
       if (me == 0) then
       print *,' in compns_stochy'
-      print*,'skeb=',skeb
       endif
 
 ! PJP stochastic physics additions
@@ -326,4 +326,156 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
       return
       end subroutine compns_stochy
 
+      subroutine compns_stochy_ocn (deltim,iret)
+!$$$  Subprogram Documentation Block
+!
+! Subprogram:  compns     Check and compute namelist frequencies
+!   Prgmmr: Iredell       Org: NP23          Date: 1999-01-26
+!
+! Abstract: This subprogram checks global spectral model namelist
+!           frequencies in hour units for validity.  If they are valid,
+!           then the frequencies are computed in timestep units.
+!           The following rules are applied:
+!             1. the timestep must be positive;
+!
+! Program History Log:
+!   2016-10-11  Phil Pegion  make the stochastic physics stand alone
+!
+! Usage:    call compns_stochy (me,deltim,nlunit, stochy_namelist,iret)
+!   Input Arguments:
+!     deltim   - real timestep in seconds
+!   Output Arguments:
+!     iret     - integer return code (0 if successful or
+!                between 1 and 8 for which rule above was broken)
+!     stochy_namelist
+!
+! Attributes:
+!   Language: Fortran 90
+!
+!$$$
+
+
+      use stochy_namelist_def
+      use mpp_mod ,only: mpp_pe,mpp_root_pe
+
+      implicit none
+
+
+      real,                 intent(in)  :: deltim
+      integer,              intent(out) :: iret
+      real tol,l_min
+      real :: rerth,circ
+      integer k,ios,nlunit
+      integer,parameter :: four=4
+
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+!
+      namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, &
+      iseed_shum,iseed_sppt,shum,shum_tau, &
+      shum_lscale,stochini,skeb_varspect_opt,sppt_sfclimit, &
+      skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, &
+      skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,&
+      shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, &
+      epbl,epbl_lscale,epbl_tau,iseed_epbl,                                    &
+      ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt
+
+      namelist /nam_sfcperts/lndp_type,lndp_var_list, lndp_prt_list, iseed_lndp, & 
+      lndp_tau,lndp_lscale 
+
+
+      rerth  =6.3712e+6      ! radius of earth (m)
+      tol=0.01  ! tolerance for calculations
+      nlunit=322
+!     spectral resolution defintion
+      ntrunc=-999
+      lon_s=-999
+      lat_s=-999
+      ! can specify up to 5 values for the stochastic physics parameters
+      ! (each is an array of length 5)
+      epbl             = -999.  ! stochastic physics tendency amplitude
+      ocnsppt          = -999.  ! stochastic physics tendency amplitude
+! logicals
+      pert_epbl = .false.
+      do_ocnsppt = .false.
+      new_lscale = .false.
+      epblint          = 0
+      ocnspptint       = 0
+      epbl_tau         = -999.  ! time scales
+      ocnsppt_tau      = -999.  ! time scales
+      epbl_lscale      = -999.  ! length scales
+      ocnsppt_lscale   = -999.  ! length scales
+      iseed_epbl       = 0      ! random seeds (if 0 use system clock)
+      iseed_epbl2      = 0      ! random seeds (if 0 use system clock)
+      iseed_ocnsppt    = 0      ! random seeds (if 0 use system clock)
+      rewind (nlunit)
+      open (unit=nlunit, file='input.nml', action='READ', status='OLD', iostat=ios)
+      read(nlunit,nam_stochy)
+
+      if (mpp_pe()==mpp_root_pe()) then
+      print *,' in compns_stochy_ocn'
+      endif
+
+! PJP stochastic physics additions
+      IF (epbl(1) > 0 ) THEN
+        pert_epbl=.true.
+      ENDIF
+      IF (ocnsppt(1) > 0 ) THEN
+        do_ocnsppt=.true.
+      ENDIF
+!    compute frequencty to update random pattern
+      IF (epblint == 0.) epblint=deltim
+      nsepbl=nint(epblint/deltim)                              ! epblint in seconds
+      IF(nsepbl<=0 .or. abs(nsepbl-epblint/deltim)>tol) THEN
+         WRITE(0,*) "ePBL interval is invalid",epblint
+        iret=9
+        return
+      ENDIF
+      IF (ocnspptint == 0.) ocnspptint=deltim
+      nsocnsppt=nint(ocnspptint/deltim)                         ! ocnspptint in seconds
+      IF(nsocnsppt<=0 .or. abs(nsocnsppt-ocnspptint/deltim)>tol) THEN
+         WRITE(0,*) "ePBL interval is invalid",ocnspptint
+        iret=9
+        return
+      ENDIF
+!calculate ntrunc if not supplied
+     if (ntrunc .LT. 1) then  
+        if (mpp_pe()==mpp_root_pe()) print*,'ntrunc not supplied, calculating'
+        circ=2*3.1415928*rerth ! start with lengthscale that is circumference of the earth
+        l_min=circ
+        do k=1,5
+           if (epbl(k).GT.0) l_min=min(epbl_lscale(k),l_min)
+           if (ocnsppt(k).GT.0) l_min=min(ocnsppt_lscale(k),l_min)
+       enddo
+       !ntrunc=1.5*circ/l_min
+       ntrunc=circ/l_min
+       if (mpp_pe()==mpp_root_pe()) print*,'ntrunc calculated from l_min',l_min,ntrunc
+     endif
+     ! ensure lat_s is a mutiple of 4 with a reminader of two
+     ntrunc=INT((ntrunc+1)/four)*four+2
+     if (mpp_pe()==mpp_root_pe()) print*,'NOTE ntrunc adjusted for even nlats',ntrunc
+
+! set up gaussian grid for ntrunc if not already defined. 
+     if (lon_s.LT.1 .OR. lat_s.LT.1) then
+        lat_s=ntrunc*1.5+1
+        lon_s=lat_s*2+4
+! Grid needs to be larger since interpolation is bi-linear
+        lat_s=lat_s*2
+        lon_s=lon_s*2
+        if (mpp_pe()==mpp_root_pe()) print*,'gaussian grid not set, defining here',lon_s,lat_s
+     endif
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+!
+!  All checks are successful.
+!
+      if (mpp_pe()==mpp_root_pe()) then
+         print *, 'ocean stochastic physics'
+         print *, ' pert_epbl : ', pert_epbl
+         print *, ' do_ocnsppt : ', do_ocnsppt
+      endif
+      iret = 0
+      if (iseed_epbl(1) > 0) iseed_epbl2(1)=iseed_epbl(1)-1234567
+!
+      return
+      end subroutine compns_stochy_ocn
+
 end module compns_stochy_mod

dezouv_stochy.f

@@ -14,7 +14,6 @@ subroutine dezouv_stochy(dev,zod,uev,vod,epsedn,epsodn,
 cc
 
 cc
-      use stochy_resol_def
       use spectral_layout_mod
       use kinddef
       implicit none

dozeuv_stochy.f

@@ -12,7 +12,6 @@ module dozeuv_stochy_mod
       subroutine dozeuv_stochy(dod,zev,uod,vev,epsedn,epsodn,
      &                  snnp1ev,snnp1od,ls_node)
 cc
-      use stochy_resol_def
       use spectral_layout_mod
       use kinddef
       implicit none

epslon_stochy.f

@@ -10,7 +10,6 @@ module epslon_stochy_mod
       subroutine epslon_stochy(epse,epso,epsedn,epsodn,
      &                  ls_node)
 cc
-      use stochy_resol_def
       use spectral_layout_mod
       use kinddef
       implicit none