commits to Psd/develop

README.standalone

@@ -0,0 +1,33 @@
+Instructions on how to build and run the stand alone ca_global code on hera
+1.  Copy enitre directory strucure
+  cd /scratch2/BMC/gsienkf/Philip.Pegion
+  cp -r NEMSfv3gfs_develop  <your path>
+
+2. Optional Compile full model: (I alredy built is, so you should not have to)
+    cd <your path>/tests/
+    ./compile.sh $PWD/../FV3 hera.intel '32BIT=Y HYDRO=N' 1 NO NO
+
+4.  compile standalone model  
+    cd ../stochastic_physics
+   # my loaded modules:
+    1) intel/18.0.5.274   3) netcdf/4.7.0   5) contrib       7) anaconda/latest
+    2) impi/2018.0.4      4) hdf5/1.10.5    6) grads/2.2.1
+
+
+    compile_standalone_ca OR compile_standalone
+
+5.  There is an input.nml which is currently set to run at C96, and the INPUT directory points to
+    C96 grid files. If you want to run at C384, you will need to rm INPUT and link in the C384 
+    directory (ln -sf INPUT_C384 INPUT) and make npx and npy 385 in input.nml
+
+5.  get an interactive session
+
+    salloc --x11=first -q debug -t 0:30:00 --nodes=1 -A <your account> -n 6
+
+6.  run job
+    srun standalone_ca OR standalone_stochy
+
+7.  post-process job to 360x180 grid (script is on hera in ~Philip.Pegion)
+    run_fregrid.sh (OR run_fregrid_c384.sh)
+
+    output file is ca_out.nc

compile_standalone

@@ -38,6 +38,4 @@ if [ $compile_all -eq 1 ];then
    $FC ${FLAGS64} stochastic_physics.F90
    ar rv libstochastic_physics.a *.o
 fi
-#$FC ${FLAGS64} get_stochy_pattern.F90
-#ar rv libstochastic_physics.a *.o
 $FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf

compile_standalone_ca

@@ -0,0 +1,17 @@
+#!/bin/sh
+compile_all=1
+if [ $compile_all -eq 1 ];then
+   rm -f *.i90 *.i *.o *.mod lib*a
+   mpif90 -C -traceback -real-size 64 -c ../FV3/gfsphysics/physics/machine.F
+   mpif90 -C -traceback -real-size 64 -c ../FV3/gfsphysics/physics/mersenne_twister.f
+   mpif90 -C -traceback -I../FMS/include -c ../FMS/platform/platform.F90
+   mpif90 -C -traceback -nowarn -DGFS_PHYS -I../FMS/include -c ../FMS/constants/constants.F90
+   mpif90 -C -traceback -I../FV3/atmos_cubed_sphere -I../FMS/include -I../FMS/fv3gfs -c fv_control_stub.F90
+   mpif90 -C -traceback -I../FV3/atmos_cubed_sphere -I../FMS/include -I../FMS/fv3gfs -c atmosphere_stub.F90
+   mpif90 -C -traceback -c -real-size 64 standalone_stochy_module.F90
+   mpif90 -C -traceback -I. -real-size 64 -c plumes.F90
+   mpif90 -DSTOCHY_UNIT_TEST -real-size 64 -C -traceback -I../FMS/fv3gfs -I../FMS/FMS_INSTALL -I../FV3/atmos_cubed_sphere -c update_ca.F90
+   mpif90 -DSTOCHY_UNIT_TEST -real-size 64 -C -traceback -I../FMS/fv3gfs -I../FMS/FMS_INSTALL -I../FV3/atmos_cubed_sphere -c cellular_automata_global.F90
+   ar rv libcellular_automata.a *.o
+fi
+mpif90 -traceback -real-size 64 -qopenmp -o standalone_ca standalone_ca.F90 -I../FV3/atmos_cubed_sphere -I../FMS/FMS_INSTALL -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lcellular_automata -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf

compns_stochy.F90

@@ -44,8 +44,10 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
       character(len=*),     intent(in)  :: input_nml_file(sz_nml)
       character(len=64),    intent(in)  :: fn_nml
       real,                 intent(in)  :: deltim
-      real tol
+      real tol,l_min
+      real :: rerth,circ,tmp_lat
       integer k,ios
+      integer,parameter :: four=4
 
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -58,6 +60,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
       namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts
       pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land
 
+      rerth  =6.3712e+6      ! radius of earth (m)
       tol=0.01  ! tolerance for calculations
 !     spectral resolution defintion
       ntrunc=-999
@@ -120,7 +123,7 @@ 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 hour to dump random patterns
+      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
@@ -207,6 +210,33 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
           pertlai(1) > 0 .OR. pertvegf(1) > 0 .OR. pertalb(1) > 0) THEN
         do_sfcperts=.true.
       ENDIF
+!calculate ntrunc if not supplied
+     if (ntrunc .LT. 1) then  
+        if (me==0) 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 (sppt(k).GT.0) l_min=min(sppt_lscale(k),l_min)
+           if (shum(k).GT.0) l_min=min(shum_lscale(k),l_min)
+           if (skeb(k).GT.0) l_min=min(skeb_lscale(k),l_min)
+       enddo
+       !ntrunc=1.5*circ/l_min
+       ntrunc=circ/l_min
+       if (me==0) 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 (me==0) 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 (me==0) print*,'gaussian grid not set, defining here',lon_s,lat_s
+     endif
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
 !  All checks are successful.

get_stochy_pattern.F90

@@ -81,6 +81,7 @@ subroutine get_random_pattern_fv3(rpattern,npatterns,&
   enddo
 
    call mp_reduce_sum(workg,lonf,latg)
+
 
 ! interpolate to cube grid
    allocate(rslmsk(lonf,latg))

standalone_ca.F90

@@ -0,0 +1,265 @@
+program  standalone_stochy_new
+
+use standalone_stochy_module
+
+use fv_mp_mod,           only: mp_start, domain_decomp
+use atmosphere_stub_mod, only: Atm,atmosphere_init_stub
+use fv_arrays_mod,       only: fv_atmos_type
+use fv_control_mod,      only: setup_pointers
+!use mpp_domains
+use mpp_mod,             only: mpp_set_current_pelist,mpp_get_current_pelist,mpp_init,mpp_pe,mpp_npes ,mpp_declare_pelist
+use mpp_domains_mod,     only: mpp_broadcast_domain,MPP_DOMAIN_TIME,mpp_domains_init ,mpp_domains_set_stack_size
+use fms_mod,             only:  fms_init
+!use time_manager_mod,    only: time_type
+use xgrid_mod,           only: grid_box_type
+use netcdf
+
+
+implicit none
+type(GFS_control_type)  :: Model
+integer, parameter      :: nlevs=64
+integer                 :: ntasks,fid,ct
+integer                 :: nthreads,omp_get_num_threads
+integer                 :: ncid,xt_dim_id,yt_dim_id,time_dim_id,xt_var_id,yt_var_id,time_var_id,ca_out_id
+integer                 :: ca1_id,ca2_id,ca3_id
+character*1             :: strid
+type(GFS_grid_type),allocatable     :: Grid(:)
+type(GFS_diag_type),allocatable :: Diag(:)
+type(GFS_statein_type),allocatable :: Statein(:)
+type(GFS_coupling_type),allocatable :: Coupling(:)
+include 'mpif.h'
+include 'netcdf.inc'
+real(kind=4) :: ts,undef
+
+integer     :: cres,blksz,nblks,ierr,my_id,i,j,nx2,ny2,nx,ny,id
+integer,target :: npx,npy
+integer     :: ng,layout(2),io_layout(2),commID,grid_type,ntiles
+integer :: halo_update_type = 1
+logical,target :: nested
+integer  :: pe,npes,stackmax=4000000
+
+real(kind=4),allocatable,dimension(:,:) :: workg
+real(kind=4),allocatable,dimension(:) :: grid_xt,grid_yt
+real(kind=8),pointer    ,dimension(:,:) :: area
+type(grid_box_type)           :: grid_box
+!type(time_type)               :: Time               ! current time
+!type(time_type)               :: Time_step          ! atmospheric time step.
+!type(time_type)               :: Time_init          ! reference time.
+!---cellular automata control parameters
+integer              :: nca             !< number of independent cellular automata
+integer              :: nlives          !< cellular automata lifetime
+integer              :: ncells          !< cellular automata finer grid
+real                 :: nfracseed       !< cellular automata seed probability
+integer              :: nseed           !< cellular automata seed frequency
+logical              :: do_ca           !< cellular automata main switch
+logical              :: ca_sgs          !< switch for sgs ca
+logical              :: ca_global       !< switch for global ca
+logical              :: ca_smooth       !< switch for gaussian spatial filter
+logical              :: isppt_deep      !< switch for combination with isppt_deep. OBS! Switches off SPPT on other tendencies!
+logical              :: isppt_pbl
+logical              :: isppt_shal
+logical              :: pert_flux   
+logical              :: pert_trigger
+integer              :: iseed_ca        !< seed for random number generation in ca scheme
+integer              :: nspinup         !< number of iterations to spin up the ca
+integer              :: ca_amplitude
+real                 :: nthresh         !< threshold used for perturbed vertical velocity
+
+NAMELIST /gfs_physics_nml/ nca, ncells, nlives, nfracseed,nseed, nthresh, &
+         do_ca,ca_sgs, ca_global,iseed_ca,ca_smooth,isppt_pbl,isppt_shal,isppt_deep,nspinup,&
+         pert_trigger,pert_flux,ca_amplitude
+
+! default values
+nca            = 1
+ncells         = 5
+nlives         = 10
+nfracseed      = 0.5
+nseed          = 100000
+iseed_ca       = 0
+nspinup        = 1
+do_ca          = .false.
+ca_sgs         = .false.
+ca_global      = .false.
+ca_smooth      = .false.
+isppt_deep     = .false.
+isppt_shal     = .false.
+isppt_pbl      = .false.
+pert_trigger   = .false.
+pert_flux      = .false.
+ca_amplitude   = 500.
+nthresh        = 0.0
+
+! open namelist file
+open (unit=565, file='input.nml', READONLY, status='OLD', iostat=ierr)
+read(565,gfs_physics_nml)
+close(565)
+! define stuff
+ng=3  ! ghost region
+undef=9.99e+20
+
+! initialize fms
+call fms_init()
+call mpp_init()
+call fms_init
+my_id=mpp_pe()
+
+call atmosphere_init_stub (grid_box, area)
+!define domain
+isd=Atm(1)%bd%isd
+ied=Atm(1)%bd%ied
+jsd=Atm(1)%bd%jsd
+jed=Atm(1)%bd%jed
+isc=Atm(1)%bd%isc
+iec=Atm(1)%bd%iec
+jsc=Atm(1)%bd%jsc
+jec=Atm(1)%bd%jec
+nx=Atm(1)%npx-1
+ny=Atm(1)%npy-1
+allocate(workg(nx,ny))
+
+! for this simple test, nblocks = ny, blksz=ny
+nblks=ny
+blksz=nx
+nthreads = omp_get_num_threads()
+Model%me=my_id
+
+Model%nca            = nca
+Model%ncells         = ncells
+Model%nlives         = nlives
+Model%nfracseed      = nfracseed
+Model%nseed          = nseed  
+Model%iseed_ca       = iseed_ca
+Model%nspinup        = nspinup
+Model%do_ca          = do_ca
+Model%ca_sgs         = ca_sgs  
+Model%ca_global      = ca_global
+Model%ca_smooth      = ca_smooth
+Model%isppt_deep     = isppt_deep
+Model%isppt_pbl      = isppt_pbl 
+Model%isppt_shal     = isppt_shal
+Model%pert_flux      = pert_flux 
+Model%pert_trigger   = pert_trigger
+Model%nthresh        = nthresh
+
+! setup GFS_init parameters
+
+!define model grid
+
+allocate(grid_xt(nx),grid_yt(ny))
+do i=1,nx
+  grid_xt(i)=i
+enddo
+do i=1,ny
+  grid_yt(i)=i
+enddo
+
+!setup GFS_coupling
+allocate(Diag(nblks))
+allocate(Coupling(nblks))
+allocate(Statein(nblks))
+write(strid,'(I1.1)') my_id+1
+fid=30+my_id
+ierr=nf90_create('ca_out.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid)
+ierr=NF90_DEF_DIM(ncid,"grid_xt",nx,xt_dim_id)
+ierr=NF90_DEF_DIM(ncid,"grid_yt",ny,yt_dim_id)
+ierr=NF90_DEF_DIM(ncid,"time",NF90_UNLIMITED,time_dim_id)
+  !> - Define the dimension variables.
+ierr=NF90_DEF_VAR(ncid,"grid_xt",NF90_FLOAT,(/ xt_dim_id /), xt_var_id)
+ierr=NF90_PUT_ATT(ncid,xt_var_id,"long_name","T-cell longitude")
+ierr=NF90_PUT_ATT(ncid,xt_var_id,"cartesian_axis","X")
+ierr=NF90_PUT_ATT(ncid,xt_var_id,"units","degrees_E")
+ierr=NF90_DEF_VAR(ncid,"grid_yt",NF90_FLOAT,(/ yt_dim_id /), yt_var_id)
+ierr=NF90_PUT_ATT(ncid,yt_var_id,"long_name","T-cell latitude")
+ierr=NF90_PUT_ATT(ncid,yt_var_id,"cartesian_axis","Y")
+ierr=NF90_PUT_ATT(ncid,yt_var_id,"units","degrees_N")
+ierr=NF90_DEF_VAR(ncid,"time",NF90_FLOAT,(/ time_dim_id /), time_var_id)
+ierr=NF90_PUT_ATT(ncid,time_var_id,"long_name","time")
+ierr=NF90_PUT_ATT(ncid,time_var_id,"units","hours since 2014-08-01 00:00:00")
+ierr=NF90_PUT_ATT(ncid,time_var_id,"cartesian_axis","T")
+ierr=NF90_PUT_ATT(ncid,time_var_id,"calendar_type","JULIAN")
+ierr=NF90_PUT_ATT(ncid,time_var_id,"calendar","JULIAN")
+!ierr=NF90_DEF_VAR(ncid,"ca_out",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca_out_id)
+!ierr=NF90_PUT_ATT(ncid,ca_out_id,"long_name","random pattern")
+!ierr=NF90_PUT_ATT(ncid,ca_out_id,"units","None")
+!ierr=NF90_PUT_ATT(ncid,ca_out_id,"missing_value",undef)
+!ierr=NF90_PUT_ATT(ncid,ca_out_id,"_FillValue",undef)
+!ierr=NF90_PUT_ATT(ncid,ca_out_id,"cell_methods","time: point")
+ierr=NF90_DEF_VAR(ncid,"ca1",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca1_id)
+ierr=NF90_PUT_ATT(ncid,ca1_id,"long_name","random pattern")
+ierr=NF90_PUT_ATT(ncid,ca1_id,"units","None")
+ierr=NF90_PUT_ATT(ncid,ca1_id,"missing_value",undef)
+ierr=NF90_PUT_ATT(ncid,ca1_id,"_FillValue",undef)
+ierr=NF90_PUT_ATT(ncid,ca1_id,"cell_methods","time: point")
+ierr=NF90_DEF_VAR(ncid,"ca2",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca2_id)
+ierr=NF90_PUT_ATT(ncid,ca2_id,"long_name","random pattern")
+ierr=NF90_PUT_ATT(ncid,ca2_id,"units","None")
+ierr=NF90_PUT_ATT(ncid,ca2_id,"missing_value",undef)
+ierr=NF90_PUT_ATT(ncid,ca2_id,"_FillValue",undef)
+ierr=NF90_PUT_ATT(ncid,ca2_id,"cell_methods","time: point")
+ierr=NF90_DEF_VAR(ncid,"ca3",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca3_id)
+ierr=NF90_PUT_ATT(ncid,ca3_id,"long_name","random pattern")
+ierr=NF90_PUT_ATT(ncid,ca3_id,"units","None")
+ierr=NF90_PUT_ATT(ncid,ca3_id,"missing_value",undef)
+ierr=NF90_PUT_ATT(ncid,ca3_id,"_FillValue",undef)
+ierr=NF90_PUT_ATT(ncid,ca3_id,"cell_methods","time: point")
+ierr=NF90_ENDDEF(ncid)
+ierr=NF90_PUT_VAR(ncid,xt_var_id,grid_xt)
+ierr=NF90_PUT_VAR(ncid,yt_var_id,grid_yt)
+! allocate diagnostics
+DO i =1,nblks
+   allocate(Diag(i)%ca_out(blksz))
+   allocate(Diag(i)%ca_deep(blksz))
+   allocate(Diag(i)%ca_turb(blksz))
+   allocate(Diag(i)%ca_shal(blksz))
+   allocate(Diag(i)%ca_rad(blksz))
+   allocate(Diag(i)%ca_micro(blksz))
+   allocate(Diag(i)%ca1(blksz))
+   allocate(Diag(i)%ca2(blksz))
+   allocate(Diag(i)%ca3(blksz))
+! allocate coupling
+   allocate(Coupling(i)%cape(blksz))
+   allocate(Coupling(i)%ca_out(blksz))
+   allocate(Coupling(i)%ca_deep(blksz))
+   allocate(Coupling(i)%ca_turb(blksz))
+   allocate(Coupling(i)%ca_shal(blksz))
+   allocate(Coupling(i)%ca_rad(blksz))
+   allocate(Coupling(i)%ca_micro(blksz))
+   allocate(Coupling(i)%ca1(blksz))
+   allocate(Coupling(i)%ca2(blksz))
+   allocate(Coupling(i)%ca3(blksz))
+! allocate coupling
+   allocate(Statein(i)%pgr(blksz))
+   allocate(Statein(i)%qgrs(blksz,nlevs,1))
+   allocate(Statein(i)%vvl(blksz,nlevs))
+   allocate(Statein(i)%prsl(blksz,nlevs))
+ENDDO
+ct=1
+do i=1,600
+   ts=i/8.0  ! hard coded to write out hourly based on a 450 second time-step
+   call cellular_automata_global(i-1, Statein, Coupling, Diag, &
+                          nblks, Model%levs, Model%nca, Model%ncells,          &
+                          Model%nlives, Model%nfracseed, Model%nseed,                    &
+                          Model%nthresh, Model%ca_global, Model%ca_sgs,                  &
+                          Model%iseed_ca, Model%ca_smooth, Model%nspinup,                &
+                          blksz)
+   if (mod(i,8).EQ.0) then
+      do j=1,ny
+         workg(:,j)=Diag(j)%ca1(:)   
+      enddo
+      ierr=NF90_PUT_VAR(ncid,ca1_id,workg,(/1,1,ct/))
+      do j=1,ny
+         workg(:,j)=Diag(j)%ca2(:)   
+      enddo
+      ierr=NF90_PUT_VAR(ncid,ca2_id,workg,(/1,1,ct/))
+      do j=1,ny
+         workg(:,j)=Diag(j)%ca3(:)   
+      enddo
+      ierr=NF90_PUT_VAR(ncid,ca3_id,workg,(/1,1,ct/))
+      ierr=NF90_PUT_VAR(ncid,time_var_id,ts,(/ct/))
+      ct=ct+1
+   if (my_id.EQ.0) write(6,fmt='(a,i5,4f6.3)') 'ca=',i,Diag(1)%ca1(1:4)
+   endif
+enddo
+!close(fid)
+ierr=NF90_CLOSE(ncid)
+end