diff --git a/CMakeLists.txt b/CMakeLists.txt index 3ad04d6a..3d94ec2d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,6 +53,7 @@ add_library( ./cellular_automata_global.F90 ./cellular_automata_sgs.F90 ./update_ca.F90 + ./lndp_apply_perts.F90 ) target_link_libraries(stochastic_physics sp::sp_d) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 4914c468..b61e8945 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -61,8 +61,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) 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 - namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts - pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land + 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 @@ -75,27 +75,30 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) sppt = -999. ! stochastic physics tendency amplitude shum = -999. ! stochastic boundary layer spf hum amp skeb = -999. ! stochastic KE backscatter amplitude - ! mg, sfcperts - pertz0 = -999. ! momentum roughness length amplitude - pertshc = -999. ! soil hydraulic conductivity amp - pertzt = -999. ! mom/heat roughness length amplitude - pertlai = -999. ! leaf area index amplitude - pertvegf = -999. ! vegetation fraction amplitude - pertalb = -999. ! albedo perturbations amplitude + lndp_var_list = 'XXX' + lndp_prt_list = -999. ! logicals do_sppt = .false. use_zmtnblck = .false. new_lscale = .false. do_shum = .false. do_skeb = .false. - ! mg, sfcperts - do_sfcperts = .false. - sppt_land = .false. - nsfcpert = 0 -! for sfcperts random patterns - sfc_lscale = -999. ! length scales - sfc_tau = -999. ! time scales - iseed_sfc = 0 ! random seeds (if 0 use system clock) +! C. Draper July 2020. +! input land pert variables: +! LNDP_TYPE = 0 +! no explicit land perturbations +! LNDP_Type = 1 +! this is the initial land sfc pert scheme, introduced and tested for impact on GEFS forecasts. +! see https://journals.ametsoc.org/doi/full/10.1175/MWR-D-18-0057.1 +! perturbations are assigned once at the start of the forecast +! LNDP_TYPE = 2 +! this is the newer land pert scheme, introduced and tested for impact on UFS/GDAS cycling stsyem +! perturbations are assigned at each time step (for state variables), or each time parameters are updated +! and the perturbations evolve over time. + lndp_type = 0 ! + lndp_lscale = -999. ! length scales + lndp_tau = -999. ! time scales + iseed_lndp = 0 ! random seeds (if 0 use system clock) ! for SKEB random patterns. skeb_vfilt = 0 skebint = 0 @@ -209,11 +212,6 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) iret=9 return ENDIF -! mg, sfcperts - IF (pertz0(1) > 0 .OR. pertshc(1) > 0 .OR. pertzt(1) > 0 .OR. & - 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' @@ -224,6 +222,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) 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 + if (lndp_type.GT.0) l_min=min(lndp_lscale(1),l_min) !ntrunc=1.5*circ/l_min ntrunc=circ/l_min if (me==0) print*,'ntrunc calculated from l_min',l_min,ntrunc @@ -241,7 +240,76 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) lon_s=lon_s*2 if (me==0) print*,'gaussian grid not set, defining here',lon_s,lat_s endif + +! +! land perts - parse nml input ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + select case (lndp_type) + case (0) + if (me==0) print*, & + 'no land perturbations selected' + case (1,2) + ! count requested pert variables + n_var_lndp= 0 + do k =1,size(lndp_var_list) + if ( (lndp_var_list(k) .EQ. 'XXX') .or. (lndp_prt_list(k) .LE. 0.) ) then + cycle + else + n_var_lndp=n_var_lndp+1 + lndp_var_list( n_var_lndp) = lndp_var_list(k) ! for lndp_type==2: + ! for state variables, unit is pert per hour + ! for parmaters, no time dimension in unit + ! since perturbations do not accumulate + ! (i.e., global_cycle overwrites the paramaters + ! each time it's called, so any previous perturbations + ! are lost). + lndp_prt_list( n_var_lndp) = lndp_prt_list(k) + endif + enddo + if (n_var_lndp > max_n_var_lndp) then + print*, 'ERROR: land perturbation requested for too many parameters', & + 'increase max_n_var_lndp' + iret = 10 + return + endif + + + if (lndp_type==1) then + if (me==0) print*, & + 'lndp_type=1, land perturbations will be applied to selected paramaters, using older scheme designed for S2S fcst spread' + ! sanity-check requested input + do k =1,n_var_lndp + select case (lndp_var_list(k)) + case('rz0','rzt','shc','lai','vgf','alb') + if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k) + case default + print*, 'ERROR: land perturbation requested for unknown parameter', lndp_var_list(k) + iret = 10 + return + end select + enddo + elseif(lndp_type==2) then + if (me==0) print*, & + 'land perturbations will be applied to selected paramaters, using newer scheme designed for DA ens spread' + do k =1,n_var_lndp + select case (lndp_var_list(k)) + case('vgf','smc','stc') + if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k) + case default + print*, 'ERROR: land perturbation requested for new parameter - will need to be coded in lndp_apply_pert', lndp_var_list(k) + iret = 10 + return + end select + enddo + endif + + case default + if (me==0) print*, & + 'lndp_type out of range, set to 0 (none), 1 (for fcst spread), 2 (for cycling DA spread)' + iret = 10 + return + end select ! ! All checks are successful. ! @@ -250,7 +318,8 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) print *, ' do_sppt : ', do_sppt print *, ' do_shum : ', do_shum print *, ' do_skeb : ', do_skeb - print *, ' do_sfcperts : ', do_sfcperts + print *, ' lndp_type : ', lndp_type + if (lndp_type .NE. 0) print *, ' n_var_lndp : ', n_var_lndp endif iret = 0 ! diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index 2e453f34..b3f712fd 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -6,11 +6,11 @@ module get_stochy_pattern_mod use spectral_layout_mod, only : ipt_lats_node_a, lat1s_a, lats_dim_a, & lats_node_a, lon_dim_a, len_trie_ls, & len_trio_ls, ls_dim, nodes, stochy_la2ga - use stochy_namelist_def, only : nsfcpert, ntrunc, stochini + use stochy_namelist_def, only : n_var_lndp, ntrunc, stochini use stochy_data_mod, only : gg_lats, gg_lons, inttyp, nskeb, nshum, nsppt, & rad2deg, rnlat, rpattern_sfc, rpattern_skeb, & rpattern_shum, rpattern_sppt, skebu_save, & - skebv_save, skeb_vwts, skeb_vpts, wlon + skebv_save, skeb_vwts, skeb_vpts, wlon, nlndp use stochy_gg_def, only : coslat_a use stochy_patterngenerator_mod, only: random_pattern, ndimspec, & patterngenerator_advance @@ -25,7 +25,7 @@ module get_stochy_pattern_mod private public get_random_pattern_fv3,get_random_pattern_fv3_vect - public get_random_pattern_sfc_fv3 + public get_random_pattern_fv3_sfc public dump_patterns logical :: first_call=.true. contains @@ -102,16 +102,17 @@ end subroutine get_random_pattern_fv3 !>@brief The subroutine 'get_random_pattern_fv3_sfc' converts spherical harmonics to the gaussian grid then interpolates to the cubed-sphere grid once !>@details This subroutine is for a 2-D (lat-lon) scalar field -subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& +subroutine get_random_pattern_fv3_sfc(rpattern,npatterns,& gis_stochy,xlat,xlon,blksz,nblks,maxlen,pattern_3d) !\callgraph ! generate a random pattern for stochastic physics implicit none - type(random_pattern), intent(inout) :: rpattern(npatterns) - type(stochy_internal_state), target :: gis_stochy + type(random_pattern), intent(inout) :: rpattern(npatterns) + type(stochy_internal_state), target :: gis_stochy real(kind=kind_dbl_prec), intent(in) :: xlat(:,:),xlon(:,:) integer,intent(in) :: npatterns,blksz(:),nblks,maxlen + real(kind=kind_dbl_prec), intent(out) :: pattern_3d(nblks,maxlen,n_var_lndp) integer i,j,l,lat,ierr,n,nn,k,nt real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d @@ -123,16 +124,16 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a) integer kmsk0(lonf,gis_stochy%lats_node_a),len real(kind=kind_dbl_prec) :: globalvar,globalvar0 - real(kind=kind_dbl_prec) :: pattern_3d(nblks,maxlen,nsfcpert) real(kind=kind_dbl_prec) :: pattern_1d(maxlen) real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk integer :: blk - do k=1,nsfcpert + do k=1,n_var_lndp kmsk0 = 0 glolal = 0. do n=1,npatterns - if (is_master()) print *, 'Random pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k)) + call patterngenerator_advance(rpattern(n),k,.false.) + if (is_master()) print *, 'Random pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k)) call scalarspect_to_gaugrid( & rpattern(n)%spec_e(:,:,k),rpattern(n)%spec_o(:,:,k),wrk2d,& gis_stochy%ls_node,gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,& @@ -151,7 +152,7 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& enddo call mp_reduce_sum(workg,lonf,latg) - if (is_master()) print *, 'workg after mp_reduce_sum for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(workg), maxval(workg) + if (is_master()) print *, 'workg after mp_reduce_sum for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(workg), maxval(workg) ! interpolate to cube grid @@ -166,13 +167,13 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& pattern_3d(blk,:,k)=pattern_1d(:) end associate enddo - if (is_master()) print *, '3D pattern for SFC-PERTS in get_random_pattern_sfc_fv3: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k)) + if (is_master()) print *, '3D pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k)) deallocate(rslmsk) deallocate(workg) - enddo ! loop over k, nsfcpert + enddo ! loop over k, n_var_lndp -end subroutine get_random_pattern_sfc_fv3 +end subroutine get_random_pattern_fv3_sfc !>@brief The subroutine 'get_random_pattern_fv3_vect' converts spherical harmonics to a vector on gaussian grid then interpolates to the cubed-sphere grid @@ -186,6 +187,8 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& type(stochy_internal_state), target :: gis_stochy integer, intent(in) :: levs type(random_pattern), intent(inout) :: rpattern(npatterns) + real(kind=kind_dbl_prec), intent(out) :: upattern_3d(nblks,maxlen,levs) + real(kind=kind_dbl_prec), intent(out) :: vpattern_3d(nblks,maxlen,levs) real(kind=kind_evod), dimension(len_trie_ls,2,1) :: vrtspec_e,divspec_e real(kind=kind_evod), dimension(len_trio_ls,2,1) :: vrtspec_o,divspec_o @@ -193,8 +196,6 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& integer,intent(in) :: npatterns,blksz(:),nblks,maxlen integer :: blk,len - real(kind=kind_dbl_prec) :: upattern_3d(nblks,maxlen,levs) - real(kind=kind_dbl_prec) :: vpattern_3d(nblks,maxlen,levs) real(kind=kind_dbl_prec) :: pattern_1d(maxlen) real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk integer i,j,l,lat,ierr,n,nn,k,nt @@ -393,7 +394,7 @@ subroutine dump_patterns(sfile) integer :: stochlun,k,n stochlun=99 if (is_master()) then - if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0) then + if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0 .OR. nlndp > 0 ) then OPEN(stochlun,file=sfile,form='unformatted') print*,'open ',sfile,' for output' endif @@ -415,6 +416,13 @@ subroutine dump_patterns(sfile) enddo enddo endif + if (nlndp > 0) then + do n=1,nlndp + do k=1,n_var_lndp + call write_pattern(rpattern_sfc(n),k,stochlun) + enddo + enddo + endif close(stochlun) end subroutine dump_patterns !>@brief The subroutine 'write_patterns' writes out a single pattern and the seed associated with the random number sequence diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 new file mode 100644 index 00000000..dabd48ea --- /dev/null +++ b/lndp_apply_perts.F90 @@ -0,0 +1,248 @@ +module lndp_apply_perts_mod + + use kinddef, only : kind_dbl_prec + + implicit none + + private + + public :: lndp_apply_perts + + contains + +!==================================================================== +! lndp_apply_perts +!==================================================================== +! Driver for applying perturbations to sprecified land states or parameters +! Draper, July 2020. +! Note on location: requires access to namelist_soilveg + + subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & + lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmc,param_update_flag, & + smc, slc, stc, vfrac, ierr) + + implicit none + + ! intent(in) + integer, intent(in) :: blksz(:) + integer, intent(in) :: n_var_lndp, lsoil, lsm + character(len=3), intent(in) :: lndp_var_list(:) + real(kind=kind_dbl_prec), intent(in) :: lndp_prt_list(:) + real(kind=kind_dbl_prec), intent(in) :: dtf + real(kind=kind_dbl_prec), intent(in) :: sfc_wts(:,:,:) + real(kind=kind_dbl_prec), intent(in) :: xlon(:,:) + real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) + logical, intent(in) :: param_update_flag + ! true = parameters have been updated, apply perts + real(kind=kind_dbl_prec), intent(in) :: stype(:,:) + real(kind=kind_dbl_prec), intent(in) :: maxsmc(:) + + ! intent(inout) + real(kind=kind_dbl_prec), intent(inout) :: smc(:,:,:) + real(kind=kind_dbl_prec), intent(inout) :: slc(:,:,:) + real(kind=kind_dbl_prec), intent(inout) :: stc(:,:,:) + real(kind=kind_dbl_prec), intent(inout) :: vfrac(:,:) + + ! intent(out) + integer, intent(out) :: ierr + + ! local + integer :: nblks, print_i, print_nb, i, nb + integer :: this_im, v, soiltyp, k + logical :: print_flag + + real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert + + ! decrease in applied pert with depth + real(kind=kind_dbl_prec), dimension(4), parameter :: smc_vertscale = (/1.0,0.5,0.25,0.125/) + real(kind=kind_dbl_prec), dimension(4), parameter :: stc_vertscale = (/1.0,0.5,0.25,0.125/) + + ! model-dependent values, hard-wired in noah code. + real(kind=kind_dbl_prec), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/) + real(kind=kind_dbl_prec), parameter :: minsmc = 0.02 + + ierr = 0 + + if (lsm .NE. 1 ) then + write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah, ', & + ' may need to adapt variable names for a different LSM' + ierr=10 + return + endif + + nblks = size(blksz) + + call set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) + + do nb =1,nblks + do i = 1, blksz(nb) + + !if ( nint(Sfcprop(nb)%slmsk(i)) .NE. 1) cycle ! skip if not land + + !if ( ((isot == 1) .and. (soiltyp == 16)) & + ! .or.( (isot == 0) .and. (soiltyp == 9 )) ) cycle ! skip if land-ice + + if ( smc(nb,i,1) .EQ. 1.) cycle ! skip non-soil (land-ice and non-land) + ! set printing + if ( (i==print_i) .and. (nb==print_nb) ) then + print_flag = .true. + else + print_flag=.false. + endif + + do v = 1,n_var_lndp + select case (trim(lndp_var_list(v))) + !================================================================= + ! State updates - performed every cycle + !================================================================= + case('smc') + p=5. + soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc + min_bound = minsmc + max_bound = maxsmc(soiltyp) + + do k=1,lsoil + !store frozen soil moisture + tmp_sic= smc(nb,i,k) - slc(nb,i,k) + + ! perturb total soil moisture + ! factor of sldepth*1000 converts from mm to m3/m3 + pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zs_noah(k)*1000.) + pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep + ! (necessary for state vars only) + call apply_pert('smc',pert,print_flag, smc(nb,i,k),ierr,p,min_bound, max_bound) + + ! assign all of applied pert to the liquid soil moisture + slc(nb,i,k) = smc(nb,i,k) - tmp_sic + enddo + + case('stc') + + do k=1,lsoil + pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) + pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep + ! (necessary for state vars only) + call apply_pert('stc',pert,print_flag, stc(nb,i,k),ierr) + enddo + !================================================================= + ! Parameter updates - only if param_update_flag = TRUE + !================================================================= + case('vgf') ! vegetation fraction + if (param_update_flag) then + p =5. + min_bound=0. + max_bound=1. + + pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) + endif + case default + print*, & + 'ERROR: unrecognised lndp_prt_list option in lndp_apply_pert, exiting', trim(lndp_var_list(v)) + ierr = 10 + return + end select + enddo + enddo + enddo + end subroutine lndp_apply_perts + +!==================================================================== +! apply_perts +!==================================================================== +! Apply perturbations to selected land states or parameters + + subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax) + + ! intent in + logical, intent(in) :: print_flag + real(kind=kind_dbl_prec), intent(in) :: pert + character(len=*), intent(in) :: vname ! name of variable being perturbed + + real(kind=kind_dbl_prec), optional, intent(in) :: p ! flat-top paramater, 0 = no flat-top + ! flat-top function is used for bounded variables + ! to reduce the magnitude of perturbations near boundaries. + real(kind=kind_dbl_prec), optional, intent(in) :: vmin, vmax ! min,max bounds of variable being perturbed + + ! intent (inout) + real(kind=kind_dbl_prec), intent(inout) :: state + + ! intent out + integer :: ierr + + !local + real(kind=kind_dbl_prec) :: z + + if ( print_flag ) then + write(*,*) 'LNDP - applying lndp to ',vname, ', initial value', state + endif + + ! apply perturbation + if (present(p) ) then + if ( .not. (present(vmin) .and. present(vmax) )) then + ierr=20 + print*, 'error, flat-top function requires min & max to be specified' + endif + + z = -1. + 2*(state - vmin)/(vmax - vmin) ! flat-top function + state = state + pert*(1-abs(z**p)) + else + state = state + pert + endif + + if (present(vmax)) state = min( state , vmax ) + if (present(vmin)) state = max( state , vmin ) + !state = max( min( state , vmax ), vmin ) + + if ( print_flag ) then + write(*,*) 'LNDP - applying lndp to ',vname, ', final value', state + endif + + end subroutine apply_pert + + +!==================================================================== +! set_printing_nb_i +!==================================================================== +! routine to turn on print flag for selected location +! + subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) + + implicit none + + ! intent (in) + integer, intent(in) :: blksz(:) + real(kind=kind_dbl_prec), intent(in) :: xlon(:,:) + real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) + + + ! intent (out) + integer, intent(out) :: print_i, print_nb + + ! local + integer :: nblks,nb,i + real, parameter :: plon_trunc = 114.9 + real, parameter :: plat_trunc = -26.6 + real, parameter :: delta = 1. + + nblks = size(blksz) + + print_i = -9 + print_nb = -9 + do nb = 1,nblks + do i = 1,blksz(nb) + if ( (xlon(nb,i)*57.29578 > plon_trunc) .and. (xlon(nb,i)*57.29578 < plon_trunc+delta ) .and. & + (xlat(nb,i)*57.29578 > plat_trunc ) .and. (xlat(nb,i)*57.29578 < plat_trunc+delta ) ) then + print_i=i + print_nb=nb + write(*,*) 'LNDP -print flag is on', xlon(nb,i)*57.29578, xlat(nb,i)*57.29578, nb, i + return + endif + enddo + enddo + + end subroutine set_printing_nb_i + +end module lndp_apply_perts_mod + + diff --git a/makefile b/makefile index 03abd965..ba6dc07c 100644 --- a/makefile +++ b/makefile @@ -58,7 +58,8 @@ SRCS_F90 = \ ./initialize_spectral_mod.F90 \ ./cellular_automata_global.F90 \ ./cellular_automata_sgs.F90 \ - ./update_ca.F90 + ./update_ca.F90 \ + ./lndp_apply_perts.F90 SRCS_c = diff --git a/standalone_stochy.F90 b/standalone_stochy.F90 index afb8ca12..8b4b2480 100644 --- a/standalone_stochy.F90 +++ b/standalone_stochy.F90 @@ -39,7 +39,6 @@ program standalone_stochy integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb logical stochini,sppt_logit,new_lscale logical use_zmtnblck -logical sppt_land include 'mpif.h' include 'netcdf.inc' real :: ak(nlevs+1),bk(nlevs+1) @@ -71,6 +70,7 @@ program standalone_stochy integer :: nargs,ntile_out,nlunit,pe,npes,stackmax=4000000 character*80 :: fname character*1 :: ntile_out_str +integer :: iret real(kind=4),allocatable,dimension(:,:) :: workg,tile_number real(kind=4),allocatable,dimension(:,:,:) :: workg3d @@ -187,7 +187,8 @@ program standalone_stochy enddo !setup GFS_coupling allocate(Coupling(nblks)) -call init_stochastic_physics(Model, Init_parm, ntasks, nthreads) +call init_stochastic_physics(Model, Init_parm, ntasks, nthreads, iret) +if (iret .ne. 0) print *, 'ERROR init_stochastic_physics call' ! Draper - need proper error trapping here call get_outfile(fname) write(strid,'(I1.1)') my_id+1 if (ntile_out.EQ.0) write_this_tile=.true. diff --git a/standalone_stochy_module.F90 b/standalone_stochy_module.F90 index ed2038b6..39a67308 100644 --- a/standalone_stochy_module.F90 +++ b/standalone_stochy_module.F90 @@ -23,11 +23,12 @@ module standalone_stochy_module real(kind=kind_phys) :: phour !< previous forecast hour real(kind=kind_phys) :: sppt_amp !< amplitude of sppt (to go to cld scheme) integer :: kdt !< current forecast iteration - logical :: do_sppt,do_shum,do_skeb,do_sfcperts,use_zmtnblck,do_rnda - integer :: skeb_npass,nsfcpert + logical :: do_sppt,do_shum,do_skeb,use_zmtnblck,do_rnda + integer :: skeb_npass,n_var_lndp, lndp_type character(len=65) :: fn_nml !< namelist filename character(len=256),allocatable :: input_nml_file(:) !< character string containing full namelist - real(kind=kind_phys) :: pertz0(5),pertzt(5),pertshc(5),pertlai(5),pertvegf(5),pertalb(5) + real(kind=kind_phys) :: lndp_prt_list(6) ! max_n_var_lndp, max_nlndp + character(len=3) :: lndp_var_list(6) ! max_n_var_lndp !---cellular automata control parameters integer :: nca !< number of independent cellular automata integer :: nlives !< cellular automata lifetime @@ -85,7 +86,7 @@ module standalone_stochy_module real (kind=kind_phys), allocatable :: ca1 (:) real (kind=kind_phys), allocatable :: ca2 (:) real (kind=kind_phys), allocatable :: ca3 (:) - integer :: nsfcpert=6 !< number of sfc perturbations + integer :: n_var_lndp=0 !< number of land sfc variables perturbations end type GFS_coupling_type end module standalone_stochy_module diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index abf1fe72..b413a4ae 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -19,9 +19,8 @@ module stochastic_physics !allocates and polulates the necessary arrays subroutine init_stochastic_physics(levs, blksz, dtp, input_nml_file_in, fn_nml, nlunit, & - do_sppt_in, do_shum_in, do_skeb_in, do_sfcperts_in, use_zmtnblck_out, skeb_npass_out, & - nsfcpert_out, pertz0_out, pertzt_out, pertshc_out, pertlai_out, pertalb_out, pertvegf_out, & - ak, bk, nthreads, mpiroot, mpicomm) + do_sppt_in, do_shum_in, do_skeb_in, lndp_type_in, n_var_lndp_in, use_zmtnblck_out, skeb_npass_out, & + lndp_var_list_out, lndp_prt_list_out, ak, bk, nthreads, mpiroot, mpicomm, iret) !\callgraph !use stochy_internal_state_mod use stochy_data_mod, only : nshum,rpattern_shum,init_stochdata,rpattern_sppt,nsppt,rpattern_skeb,nskeb,gg_lats,gg_lons,& @@ -33,6 +32,7 @@ subroutine init_stochastic_physics(levs, blksz, dtp, input_nml_file_in, fn_nml, use mpi_wrapper, only : mpi_wrapper_initialize,mype,npes,is_master implicit none +integer, intent(out) :: iret ! Interface variables @@ -41,18 +41,18 @@ subroutine init_stochastic_physics(levs, blksz, dtp, input_nml_file_in, fn_nml, real(kind=kind_dbl_prec), intent(in) :: dtp character(len=*), intent(in) :: input_nml_file_in(:) character(len=*), intent(in) :: fn_nml -logical, intent(in) :: do_sppt_in, do_shum_in, do_skeb_in, do_sfcperts_in +logical, intent(in) :: do_sppt_in, do_shum_in, do_skeb_in +integer, intent(in) :: lndp_type_in, n_var_lndp_in +real(kind=kind_dbl_prec), intent(in) :: ak(:), bk(:) logical, intent(out) :: use_zmtnblck_out -integer, intent(out) :: skeb_npass_out, nsfcpert_out +integer, intent(out) :: skeb_npass_out +character(len=3), dimension(max_n_var_lndp), intent(out) :: lndp_var_list_out +real(kind=kind_dbl_prec), dimension(max_n_var_lndp), intent(out) :: lndp_prt_list_out -real(kind=kind_dbl_prec), intent(out) :: pertz0_out(:),pertzt_out(:),pertshc_out(:) -real(kind=kind_dbl_prec), intent(out) :: pertlai_out(:),pertalb_out(:),pertvegf_out(:) -real(kind=kind_dbl_prec), intent(in) :: ak(:), bk(:) ! Local variables real(kind=kind_dbl_prec), parameter :: con_pi =4.0d0*atan(1.0d0) integer :: nblks -integer :: iret real*8 :: PRSI(levs),PRSL(levs),dx real, allocatable :: skeb_vloc(:) integer :: k,kflip,latghf,blk,k2 @@ -75,35 +75,40 @@ subroutine init_stochastic_physics(levs, blksz, dtp, input_nml_file_in, fn_nml, gis_stochy%me=me gis_stochy%nodes=nodes call init_stochdata(levs,dtp,input_nml_file_in,fn_nml,nlunit,iret) +if (iret .ne. 0) return ! check namelist entries for consistency if (do_sppt_in.neqv.do_sppt) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_sppt and sppt' + iret = 20 return else if (do_shum_in.neqv.do_shum) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_shum and shum' + iret = 20 return else if (do_skeb_in.neqv.do_skeb) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & & ' namelist settings do_skeb and skeb' + iret = 20 + return +else if (lndp_type_in /= lndp_type) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings lndp_type in physics and nam_sfcperts' + iret = 20 return -else if (do_sfcperts_in.neqv.do_sfcperts) then ! mg, sfc-perts +else if (n_var_lndp_in /= n_var_lndp) then write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & - & ' namelist settings do_sfcperts and pertz0 / pertshc / pertzt / pertlai / pertvegf / pertalb' + & ' namelist settings n_var_lndp in physics nml, and lndp_* in nam_sfcperts' + iret = 20 return end if ! update remaining model configuration parameters from namelist use_zmtnblck_out=use_zmtnblck skeb_npass_out=skeb_npass -nsfcpert_out=nsfcpert ! mg, sfc-perts -pertz0_out=pertz0 ! mg, sfc-perts -pertzt_out=pertzt ! mg, sfc-perts -pertshc_out=pertshc ! mg, sfc-perts -pertlai_out=pertlai ! mg, sfc-perts -pertalb_out=pertalb ! mg, sfc-perts -pertvegf_out=pertvegf ! mg, sfc-perts -if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return +lndp_var_list_out=lndp_var_list +lndp_prt_list_out=lndp_prt_list +if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0) ) return allocate(sl(levs)) do k=1,levs sl(k)= 0.5*(ak(k)/101300.+bk(k)+ak(k+1)/101300.0+bk(k+1)) ! si are now sigmas @@ -214,15 +219,19 @@ end subroutine init_stochastic_physics !>@details It updates the AR(1) in spectral space !allocates and polulates the necessary arrays -subroutine run_stochastic_physics(levs, kdt, phour, blksz, xlat, xlon, sppt_wts, shum_wts, skebu_wts, skebv_wts, nthreads) +subroutine run_stochastic_physics(levs, kdt, phour, blksz, xlat, xlon, sppt_wts, shum_wts, skebu_wts, & + skebv_wts, sfc_wts,nthreads) !\callgraph !use stochy_internal_state_mod use stochy_data_mod, only : nshum,rpattern_shum,rpattern_sppt,nsppt,rpattern_skeb,nskeb,& - rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_sppt,vfact_shum,vfact_skeb -use get_stochy_pattern_mod,only : get_random_pattern_fv3,get_random_pattern_fv3_vect,dump_patterns + rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_sppt,vfact_shum,vfact_skeb, & + rpattern_sfc, nlndp +use get_stochy_pattern_mod,only : get_random_pattern_fv3,get_random_pattern_fv3_vect,dump_patterns, & + get_random_pattern_fv3_sfc use stochy_resol_def , only : latg,lonf -use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,fhstoch,nssppt,nsshum,nsskeb,sppt_logit +use stochy_namelist_def, only : do_shum,do_sppt,do_skeb,fhstoch,nssppt,nsshum,nsskeb,sppt_logit, & + lndp_type, n_var_lndp use mpi_wrapper, only: is_master use spectral_layout_mod,only:ompthreads implicit none @@ -237,137 +246,103 @@ subroutine run_stochastic_physics(levs, kdt, phour, blksz, xlat, xlon, sppt_wts, real(kind=kind_dbl_prec), intent(inout) :: shum_wts(:,:,:) real(kind=kind_dbl_prec), intent(inout) :: skebu_wts(:,:,:) real(kind=kind_dbl_prec), intent(inout) :: skebv_wts(:,:,:) +real(kind=kind_dbl_prec), intent(inout) :: sfc_wts(:,:,:) integer, intent(in) :: nthreads -real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:) +real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:), tmpl_wts(:,:,:) !D-grid integer :: k integer j,ierr,i integer :: nblks, blk, len, maxlen character*120 :: sfile character*6 :: STRFH +logical :: do_advance_pattern -if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) ) return +if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) ) return ! Update number of threads in shared variables in spectral_layout_mod and set block-related variables ompthreads = nthreads nblks = size(blksz) maxlen = maxval(blksz(:)) -! check to see if it is time to write out random patterns -if (fhstoch.GE. 0 .AND. MOD(phour,fhstoch) .EQ. 0) then - write(STRFH,FMT='(I6.6)') nint(phour) - sfile='stoch_out.F'//trim(STRFH) - call dump_patterns(sfile) -endif -allocate(tmp_wts(nblks,maxlen)) -allocate(tmpu_wts(nblks,maxlen,levs)) -allocate(tmpv_wts(nblks,maxlen,levs)) -if (do_sppt) then - if (mod(kdt,nssppt) == 1 .or. nssppt == 1) then - call get_random_pattern_fv3(rpattern_sppt,nsppt,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmp_wts) - DO blk=1,nblks - len=blksz(blk) - DO k=1,levs - sppt_wts(blk,1:len,k)=tmp_wts(blk,1:len)*vfact_sppt(k) - ENDDO - if (sppt_logit) sppt_wts(blk,:,:) = (2./(1.+exp(sppt_wts(blk,:,:))))-1. - sppt_wts(blk,:,:) = sppt_wts(blk,:,:)+1.0 - ENDDO - endif -endif -if (do_shum) then - if (mod(kdt,nsshum) == 1 .or. nsshum == 1) then - call get_random_pattern_fv3(rpattern_shum,nshum,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmp_wts) - DO blk=1,nblks - len=blksz(blk) - DO k=1,levs - shum_wts(blk,1:len,k)=tmp_wts(blk,1:len)*vfact_shum(k) - ENDDO - ENDDO - endif -endif -if (do_skeb) then - if (mod(kdt,nsskeb) == 1 .or. nsskeb == 1) then - call get_random_pattern_fv3_vect(rpattern_skeb,nskeb,gis_stochy,levs,xlat,xlon,blksz,nblks,maxlen,tmpu_wts,tmpv_wts) - DO blk=1,nblks - len=blksz(blk) - DO k=1,levs - skebu_wts(blk,1:len,k)=tmpu_wts(blk,1:len,k)*vfact_skeb(k) - skebv_wts(blk,1:len,k)=tmpv_wts(blk,1:len,k)*vfact_skeb(k) - ENDDO - ENDDO - endif -endif -deallocate(tmp_wts) -deallocate(tmpu_wts) -deallocate(tmpv_wts) + +if ( (lndp_type==1) .and. (kdt==0) ) then ! old land pert scheme called once at start + write(0,*) 'calling get_random_pattern_fv3_sfc' + allocate(tmpl_wts(nblks,maxlen,n_var_lndp)) + call get_random_pattern_fv3_sfc(rpattern_sfc,nlndp,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmpl_wts) + DO blk=1,nblks + len=blksz(blk) + ! for perturbing vars or states, saved value is N(0,1) and apply scaling later. + DO k=1,n_var_lndp + sfc_wts(blk,1:len,k) = tmpl_wts(blk,1:len,k) + ENDDO + ENDDO + deallocate(tmpl_wts) +else + ! check to see if it is time to write out random patterns + if (fhstoch.GE. 0 .AND. MOD(phour,fhstoch) .EQ. 0) then + write(STRFH,FMT='(I6.6)') nint(phour) + sfile='stoch_out.F'//trim(STRFH) + call dump_patterns(sfile) + endif + allocate(tmp_wts(nblks,maxlen)) + allocate(tmpu_wts(nblks,maxlen,levs)) + allocate(tmpv_wts(nblks,maxlen,levs)) + if (do_sppt) then + if (mod(kdt,nssppt) == 1 .or. nssppt == 1) then + call get_random_pattern_fv3(rpattern_sppt,nsppt,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmp_wts) + DO blk=1,nblks + len=blksz(blk) + DO k=1,levs + sppt_wts(blk,1:len,k)=tmp_wts(blk,1:len)*vfact_sppt(k) + ENDDO + if (sppt_logit) sppt_wts(blk,:,:) = (2./(1.+exp(sppt_wts(blk,:,:))))-1. + sppt_wts(blk,:,:) = sppt_wts(blk,:,:)+1.0 + ENDDO + endif + endif + if (do_shum) then + if (mod(kdt,nsshum) == 1 .or. nsshum == 1) then + call get_random_pattern_fv3(rpattern_shum,nshum,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmp_wts) + DO blk=1,nblks + len=blksz(blk) + DO k=1,levs + shum_wts(blk,1:len,k)=tmp_wts(blk,1:len)*vfact_shum(k) + ENDDO + ENDDO + endif + endif + if (do_skeb) then + if (mod(kdt,nsskeb) == 1 .or. nsskeb == 1) then + call get_random_pattern_fv3_vect(rpattern_skeb,nskeb,gis_stochy,levs,xlat,xlon,blksz,nblks,maxlen,tmpu_wts,tmpv_wts) + DO blk=1,nblks + len=blksz(blk) + DO k=1,levs + skebu_wts(blk,1:len,k)=tmpu_wts(blk,1:len,k)*vfact_skeb(k) + skebv_wts(blk,1:len,k)=tmpv_wts(blk,1:len,k)*vfact_skeb(k) + ENDDO + ENDDO + endif + endif + if ( lndp_type .EQ. 2 ) then + ! add time check? + allocate(tmpl_wts(nblks,maxlen,n_var_lndp)) + call get_random_pattern_fv3_sfc(rpattern_sfc,nlndp,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmpl_wts) + DO blk=1,nblks + len=blksz(blk) + ! for perturbing vars or states, saved value is N(0,1) and apply scaling later. + DO k=1,n_var_lndp + sfc_wts(blk,1:len,k) = tmpl_wts(blk,1:len,k) + ENDDO + ENDDO + deallocate(tmpl_wts) + endif + deallocate(tmp_wts) + deallocate(tmpu_wts) + deallocate(tmpv_wts) + +end if end subroutine run_stochastic_physics end module stochastic_physics - - -module stochastic_physics_sfc - -use kinddef, only : kind_dbl_prec - -implicit none - -private - -public :: run_stochastic_physics_sfc - -contains - -subroutine run_stochastic_physics_sfc(blksz, xlat, xlon, sfc_wts) - -!\callgraph -use mpi_wrapper, only : is_master -use stochy_internal_state_mod -use stochy_data_mod, only : rad2deg,INTTYP,wlon,rnlat,gis_stochy, rpattern_sfc,npsfc ! mg, sfc-perts -use get_stochy_pattern_mod,only : get_random_pattern_sfc_fv3 ! mg, sfc-perts -use stochy_resol_def , only : latg,lonf -use stochy_namelist_def, only : do_sfcperts, nsfcpert -implicit none - -! Interface variables -integer, intent(in) :: blksz(:) -real(kind=kind_dbl_prec), intent(in) :: xlat(:,:) -real(kind=kind_dbl_prec), intent(in) :: xlon(:,:) -real(kind=kind_dbl_prec), intent(out) :: sfc_wts(:,:,:) - -real,allocatable :: tmpsfc_wts(:,:,:) -!D-grid -integer :: k -integer j,ierr,i -integer :: nblks, blk, len, maxlen -character*120 :: sfile -character*6 :: STRFH - -if (.NOT. do_sfcperts) return - -! Set block-related variables -nblks = size(blksz) -maxlen = maxval(blksz(:)) - -allocate(tmpsfc_wts(nblks,maxlen,nsfcpert)) ! mg, sfc-perts -if (is_master()) then - print*,'In run_stochastic_physics_sfc' -endif -call get_random_pattern_sfc_fv3(rpattern_sfc,npsfc,gis_stochy,xlat,xlon,blksz,nblks,maxlen,tmpsfc_wts) -DO blk=1,nblks - len=blksz(blk) - DO k=1,nsfcpert - sfc_wts(blk,1:len,k)=tmpsfc_wts(blk,1:len,k) - ENDDO -ENDDO -if (is_master()) then - print*,'tmpsfc_wts(blk,1,:) =',tmpsfc_wts(1,1,1),tmpsfc_wts(1,1,2),tmpsfc_wts(1,1,3),tmpsfc_wts(1,1,4),tmpsfc_wts(1,1,5) - print*,'min(tmpsfc_wts(:,:,:)) =',minval(tmpsfc_wts(:,:,:)) -endif -deallocate(tmpsfc_wts) - -end subroutine run_stochastic_physics_sfc - -end module stochastic_physics_sfc diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index 84f46311..d835ba41 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -27,7 +27,7 @@ module stochy_data_mod integer, public :: nsppt=0 integer, public :: nshum=0 integer, public :: nskeb=0 - integer, public :: npsfc=0 + integer, public :: nlndp=0 ! this is the number of different patterns (determined by the tau/lscale input) real*8, public,allocatable :: sl(:) real(kind=kind_dbl_prec),public, allocatable :: vfact_sppt(:),vfact_shum(:),vfact_skeb(:) @@ -52,7 +52,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) character(len=64), intent(in) :: fn_nml real, intent(in) :: delt integer, intent(out) :: iret - real :: pertsfc(1) + real :: ones(5) real :: rnn1 integer :: nn,nspinup,k,nm,spinup_efolds,stochlun,ierr,n @@ -67,8 +67,9 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) iret=0 ! read in namelist call compns_stochy (me,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) + if (iret/=0) return ! need to make sure that non-zero irets are being trapped. if(is_master()) print*,'in init stochdata',nodes,lat_s - if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return + if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0) ) return ! initialize the specratl pattern generatore (including gaussian grid decomposition) ! if (nodes.GE.lat_s/2) then ! lat_s=(int(nodes/12)+1)*24 @@ -105,32 +106,23 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) endif enddo if (is_master()) print *,'nskeb = ',nskeb - ! mg, sfc-perts - do n=1,size(pertz0) - if (pertz0(n) > 0 .or. pertzt(n)>0 .or. pertshc(n)>0 .or. & - pertvegf(n)>0 .or. pertlai(n)>0 .or. pertalb(n)>0) then - npsfc=npsfc+1 - else - exit - endif - enddo - if (is_master()) then - if (npsfc > 0) then - print *,' npsfc = ', npsfc - print *,' pertz0 = ', pertz0 - print *,' pertzt = ', pertzt - print *,' pertshc = ', pertshc - print *,' pertlai = ', pertlai - print *,' pertalb = ', pertalb - print *,' pertvegf = ', pertvegf - endif - endif + ! Draper: nlndp>1 was not properly coded. Hardcode to 1 for now + !do n=1,size(lndp_z0) + ! if (lndp_z0(n) > 0 .or. lndp_zt(n)>0 .or. lndp_hc(n)>0 .or. & + ! lndp_vf(n)>0 .or. lndp_la(n)>0 .or. lndp_al(n)>0) then + ! nlndp=nlndp+1 + ! else + ! exit + ! endif + !enddo + if (n_var_lndp>0) nlndp=1 + if (is_master()) print *,' nlndp = ', nlndp if (nsppt > 0) allocate(rpattern_sppt(nsppt)) if (nshum > 0) allocate(rpattern_shum(nshum)) if (nskeb > 0) allocate(rpattern_skeb(nskeb)) ! mg, sfc perts - if (npsfc > 0) allocate(rpattern_sfc(npsfc)) + if (nlndp > 0) allocate(rpattern_sfc(nlndp)) ! if stochini is true, then read in pattern from a file if (is_master()) then @@ -148,7 +140,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) spinup_efolds = 0 if (nsppt > 0) then if (is_master()) print *, 'Initialize random pattern for SPPT' - call patterngenerator_init(sppt_lscale,spptint,sppt_tau,sppt,iseed_sppt,rpattern_sppt, & + call patterngenerator_init(sppt_lscale(1:nsppt),spptint,sppt_tau(1:nsppt),sppt(1:nsppt),iseed_sppt,rpattern_sppt, & lonf,latg,jcap,gis_stochy%ls_node,nsppt,1,0,new_lscale) do n=1,nsppt nspinup = spinup_efolds*sppt_tau(n)/spptint @@ -180,7 +172,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) endif if (nshum > 0) then if (is_master()) print *, 'Initialize random pattern for SHUM' - call patterngenerator_init(shum_lscale,shumint,shum_tau,shum,iseed_shum,rpattern_shum, & + call patterngenerator_init(shum_lscale(1:nshum),shumint,shum_tau(1:nshum),shum(1:nshum),iseed_shum,rpattern_shum, & lonf,latg,jcap,gis_stochy%ls_node,nshum,1,0,new_lscale) do n=1,nshum nspinup = spinup_efolds*shum_tau(n)/shumint @@ -215,8 +207,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) ! determine number of skeb levels to deal with temperoal/vertical correlations skeblevs=nint(skeb_tau(1)/skebint*skeb_vdof) ! backscatter noise. - if (is_master()) print *, 'Initialize random pattern for SKEB',skeblevs - call patterngenerator_init(skeb_lscale,skebint,skeb_tau,skeb,iseed_skeb,rpattern_skeb, & + call patterngenerator_init(skeb_lscale(1:nskeb),skebint,skeb_tau(1:nskeb),skeb(1:nskeb),iseed_skeb,rpattern_skeb, & lonf,latg,jcap,gis_stochy%ls_node,nskeb,skeblevs,skeb_varspect_opt,new_lscale) do n=1,nskeb do k=1,skeblevs @@ -306,38 +297,43 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) endif ! skeb > 0 ! mg, sfc-perts -if (npsfc > 0) then - pertsfc(1) = 1. - call patterngenerator_init(sfc_lscale,delt,sfc_tau,pertsfc,iseed_sfc,rpattern_sfc, & - lonf,latg,jcap,gis_stochy%ls_node,npsfc,nsfcpert,0,new_lscale) - do n=1,npsfc - if (is_master()) print *, 'Initialize random pattern for SFC-PERTS',n - do k=1,nsfcpert - nspinup = spinup_efolds*sfc_tau(n)/delt - call getnoise(rpattern_sfc(n),noise_e,noise_o) - do nn=1,len_trie_ls - rpattern_sfc(n)%spec_e(nn,1,k)=noise_e(nn,1) - rpattern_sfc(n)%spec_e(nn,2,k)=noise_e(nn,2) - nm = rpattern_sfc(n)%idx_e(nn) - if (nm .eq. 0) cycle - rpattern_sfc(n)%spec_e(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) - rpattern_sfc(n)%spec_e(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) - enddo - do nn=1,len_trio_ls - rpattern_sfc(n)%spec_o(nn,1,k)=noise_o(nn,1) - rpattern_sfc(n)%spec_o(nn,2,k)=noise_o(nn,2) - nm = rpattern_sfc(n)%idx_o(nn) - if (nm .eq. 0) cycle - rpattern_sfc(n)%spec_o(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) - rpattern_sfc(n)%spec_o(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) - enddo - do nn=1,nspinup - call patterngenerator_advance(rpattern_sfc(n),k,.false.) - enddo - if (is_master()) print *, 'Random pattern for SFC-PERTS: k, min, max ',k, minval(rpattern_sfc(1)%spec_o(:,:,k)), maxval(rpattern_sfc(1)%spec_o(:,:,k)) - enddo ! k, nsfcpert - enddo ! n, npsfc - endif ! npsfc > 0 +if (nlndp > 0) then + ones = 1. + call patterngenerator_init(lndp_lscale(1:nlndp),delt,lndp_tau(1:nlndp),ones(1:nlndp),iseed_lndp,rpattern_sfc, & + lonf,latg,jcap,gis_stochy%ls_node,nlndp,n_var_lndp,0,new_lscale) + do n=1,nlndp + if (is_master()) print *, 'Initialize random pattern for LNDP PERTS' + do k=1,n_var_lndp + nspinup = spinup_efolds*lndp_tau(n)/delt + if (stochini) then + call read_pattern(rpattern_sfc(n),k,stochlun) + if (is_master()) print *, 'lndp pattern read',n,k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k)) + else + call getnoise(rpattern_sfc(n),noise_e,noise_o) + do nn=1,len_trie_ls + rpattern_sfc(n)%spec_e(nn,1,k)=noise_e(nn,1) + rpattern_sfc(n)%spec_e(nn,2,k)=noise_e(nn,2) + nm = rpattern_sfc(n)%idx_e(nn) + if (nm .eq. 0) cycle + rpattern_sfc(n)%spec_e(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) + rpattern_sfc(n)%spec_e(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_e(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) + enddo + do nn=1,len_trio_ls + rpattern_sfc(n)%spec_o(nn,1,k)=noise_o(nn,1) + rpattern_sfc(n)%spec_o(nn,2,k)=noise_o(nn,2) + nm = rpattern_sfc(n)%idx_o(nn) + if (nm .eq. 0) cycle + rpattern_sfc(n)%spec_o(nn,1,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,1,k)*rpattern_sfc(n)%varspectrum(nm) + rpattern_sfc(n)%spec_o(nn,2,k) = rpattern_sfc(n)%stdev*rpattern_sfc(n)%spec_o(nn,2,k)*rpattern_sfc(n)%varspectrum(nm) + enddo + do nn=1,nspinup + call patterngenerator_advance(rpattern_sfc(n),k,.false.) + enddo + if (is_master()) print *, 'lndp pattern initialized, ',n, k, minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k)) + endif ! stochini + enddo ! k, n_var_lndp + enddo ! n, nlndp + endif ! nlndp > 0 if (is_master() .and. stochini) CLOSE(stochlun) deallocate(noise_e,noise_o) end subroutine init_stochdata diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index 91bb9531..310713a7 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -10,6 +10,7 @@ module stochy_namelist_def implicit none public + integer, parameter :: max_n_var_lndp = 6 ! must match value used in GFS_typedefs integer nssppt,nsshum,nsskeb,lon_s,lat_s,ntrunc ! pjp stochastic phyics @@ -29,13 +30,11 @@ module stochy_namelist_def logical use_zmtnblck logical do_shum,do_sppt,do_skeb -! mg surface perturbations - real(kind=kind_dbl_prec), dimension(5) :: sfc_lscale,sfc_tau - real(kind=kind_dbl_prec), dimension(5) :: pertz0,pertshc,pertzt - real(kind=kind_dbl_prec), dimension(5) :: pertlai,pertvegf,pertalb - integer nsfcpert - integer(8),dimension(5) ::iseed_sfc - logical sppt_land - logical do_sfcperts + real(kind=kind_dbl_prec), dimension(5) :: lndp_lscale,lndp_tau + integer n_var_lndp + integer(8),dimension(5) ::iseed_lndp + integer lndp_type + character(len=3), dimension(max_n_var_lndp) :: lndp_var_list + real(kind=kind_dbl_prec), dimension(max_n_var_lndp) :: lndp_prt_list end module stochy_namelist_def