From 1feea01e820d5b088bec3531b0ef08a46ae5a971 Mon Sep 17 00:00:00 2001 From: climbfuji Date: Tue, 9 Jul 2019 13:26:05 -0600 Subject: [PATCH] Merge portability features and bugfixes from CCPP version, make necessary changes to use stochastic physics from separate GitHub repository, trim trailing whitespaces --- cellular_automata.f90 | 102 ++--- compns_stochy.F90 | 18 +- dezouv_stochy.f | 10 +- dozeuv_stochy.f | 10 +- epslon_stochy.f | 10 +- ..._to_grid_stochy.f => four_to_grid_stochy.F | 102 ++++- get_lats_node_a_stochy.f | 14 +- get_ls_node_stochy.f | 10 +- get_stochy_pattern.F90 | 110 +++--- getcon_lag_stochy.f | 14 +- getcon_spectral.F90 | 27 +- glats_stochy.f | 7 + gozrineo_stochy.f | 10 +- initialize_spectral_mod.F90 | 64 ++-- makefile | 20 +- num_parthds_stochy.f | 2 +- pln2eo_stochy.f | 10 +- plumes.f90 | 348 +++++++++--------- setlats_a_stochy.f | 12 +- setlats_lag_stochy.f | 10 +- spectral_layout.F90 | 278 ++++++++++++++ spectral_layout.f | 28 -- stochastic_physics.F90 | 137 ++++--- stochy_data_mod.F90 | 29 +- stochy_gg_def.f | 2 +- stochy_internal_state_mod.F90 | 10 +- stochy_namelist_def.F90 | 2 +- stochy_patterngenerator.F90 | 23 +- stochy_resol_def.f | 2 +- sumfln_stochy.f | 28 +- update_ca.f90 | 108 +++--- 31 files changed, 1017 insertions(+), 540 deletions(-) rename four_to_grid_stochy.f => four_to_grid_stochy.F (75%) create mode 100644 spectral_layout.F90 delete mode 100644 spectral_layout.f diff --git a/cellular_automata.f90 b/cellular_automata.f90 index b1d98323..b8da0c29 100644 --- a/cellular_automata.f90 +++ b/cellular_automata.f90 @@ -18,20 +18,20 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !L.Bengtsson, 2017-06 !This program evolves a cellular automaton uniform over the globe given -!the flag ca_global, if instead ca_sgs is .true. it evolves a cellular automata conditioned on -!perturbed grid-box mean field. The perturbations to the mean field are given by a +!the flag ca_global, if instead ca_sgs is .true. it evolves a cellular automata conditioned on +!perturbed grid-box mean field. The perturbations to the mean field are given by a !stochastic gaussian skewed (SGS) distribution. !If ca_global is .true. it weighs the number of ca (nca) together to produce 1 output pattern -!If instead ca_sgs is given, it produces nca ca: -! 1 CA_DEEP = deep convection -! 2 CA_SHAL = shallow convection -! 3 CA_TURB = turbulence -! 4 CA_RAD = radiation -! 5 CA_MICRO = microphysics +!If instead ca_sgs is given, it produces nca ca: +! 1 CA_DEEP = deep convection +! 2 CA_SHAL = shallow convection +! 3 CA_TURB = turbulence +! 4 CA_RAD = radiation +! 5 CA_MICRO = microphysics -!PLEASE NOTE: This is considered to be version 0 of the cellular automata code for FV3GFS, some functionally -!is missing/limited. +!PLEASE NOTE: This is considered to be version 0 of the cellular automata code for FV3GFS, some functionally +!is missing/limited. integer,intent(in) :: kstep,ncells,nca,nlives,nseed,iseed_ca,nspinup real,intent(in) :: nfracseed,nthresh @@ -63,12 +63,12 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !nca :: switch for number of cellular automata to be used. !ca_global :: switch for global cellular automata -!ca_sgs :: switch for cellular automata conditioned on SGS perturbed vertvel. +!ca_sgs :: switch for cellular automata conditioned on SGS perturbed vertvel. !nfracseed :: switch for number of random cells initially seeded !nlives :: switch for maximum number of lives a cell can have !nspinup :: switch for number of itterations to spin up the ca -!ncells :: switch for higher resolution grid e.g ncells=4 -! gives 4x4 times the FV3 model grid resolution. +!ncells :: switch for higher resolution grid e.g ncells=4 +! gives 4x4 times the FV3 model grid resolution. !ca_smooth :: switch to smooth the cellular automata !nthresh :: threshold of perturbed vertical velocity used in case of sgs !nca_plumes :: compute number of CA-cells ("plumes") within a NWP gridbox. @@ -96,12 +96,12 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & stop endif - if(ca_global == .true. .and. ca_sgs == .true.)then + if(ca_global .and. ca_sgs)then write(0,*)'Namelist options ca_global and ca_sgs cannot both be true - exiting' stop endif - if(ca_sgs == .true. .and. ca_smooth == .true.)then + if(ca_sgs .and. ca_smooth)then write(0,*)'Currently ca_smooth does not work with ca_sgs - exiting' stop endif @@ -109,15 +109,15 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & call atmosphere_resolution (nlon, nlat, global=.false.) isize=nlon+2*halo jsize=nlat+2*halo - !nlon,nlat is the compute domain - without haloes - !mlon,mlat is the cubed-sphere tile size. + !nlon,nlat is the compute domain - without haloes + !mlon,mlat is the cubed-sphere tile size. inci=ncells incj=ncells - + nxc=nlon*ncells nyc=nlat*ncells - + nxch=nxc+2*halo nych=nyc+2*halo @@ -148,22 +148,22 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & allocate(CAavg(nlon,nlat)) allocate(CA_TURB(nlon,nlat)) allocate(CA_RAD(nlon,nlat)) - allocate(CA_DEEP(nlon,nlat)) + allocate(CA_DEEP(nlon,nlat)) allocate(CA_MICRO(nlon,nlat)) allocate(CA_SHAL(nlon,nlat)) allocate(noise(nxc,nyc,nca)) allocate(noise1D(nxc*nyc)) - + !Initialize: Detfield(:,:,:)=0. vertvelmean(:,:) =0. vertvelsum(:,:)=0. - cloud(:,:)=0. + cloud(:,:)=0. humidity(:,:)=0. condition(:,:)=0. cape(:,:)=0. vertvelhigh(:,:)=0. - noise(:,:,:) = 0.0 + noise(:,:,:) = 0.0 noise1D(:) = 0.0 iini(:,:,:) = 0 ilives(:,:) = 0 @@ -175,7 +175,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & CA_DEEP(:,:) = 0.0 CA_MICRO(:,:) = 0.0 CA_SHAL(:,:) = 0.0 - + !Put the blocks of model fields into a 2d array levs=nlev blocksz=blocksize @@ -187,16 +187,16 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & isc = Atm_block%isc iec = Atm_block%iec jsc = Atm_block%jsc - jec = Atm_block%jec + jec = Atm_block%jec do blk = 1,Atm_block%nblks do ix = 1, Atm_block%blksz(blk) i = Atm_block%index(blk)%ii(ix) - isc + 1 j = Atm_block%index(blk)%jj(ix) - jsc + 1 - cape(i,j) = Coupling(blk)%cape(ix) + cape(i,j) = Coupling(blk)%cape(ix) surfp(i,j) = Statein(blk)%pgr(ix) humidity(i,j)=Statein(blk)%qgrs(ix,k850,1) !about 850 hpa - do k = 1,k350 !Lower troposphere: level k350 is about 350hPa + do k = 1,k350 !Lower troposphere: level k350 is about 350hPa omega(i,j,k) = Statein(blk)%vvl(ix,k) ! layer mean vertical velocity in pa/sec pressure(i,j,k) = Statein(blk)%prsl(ix,k) ! layer mean pressure in Pa enddo @@ -235,7 +235,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & ! iseed is elapsed time since unix epoch began (secs) ! truncate to 4 byte integer count_trunc = iscale*(count/iscale) - count4 = count - count_trunc + nf*ra + count4 = count - count_trunc + nf*ra else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. @@ -254,7 +254,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !Initiate the cellular automaton with random numbers larger than nfracseed - + do j = 1,nyc do i = 1,nxc if (noise(i,j,nf) > nfracseed ) then @@ -266,7 +266,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo !nf - + !In case we want to condition the cellular automaton on a large scale field !we here set the "condition" variable to a different model field depending !on nf. (this is not used if ca_global = .true.) @@ -274,15 +274,15 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & CAstore = 0. do nf=1,nca !update each ca - - if(ca_sgs == .true.)then - if(nf==1)then + if(ca_sgs)then + + if(nf==1)then inci=ncells incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) + condition(i,j)=cape(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -317,7 +317,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & do j = 1,nyc do i = 1,nxc - ilives(i,j)=int(real(nlives)*alpha*noise(i,j,nf)) + ilives(i,j)=int(real(nlives)*alpha*noise(i,j,nf)) enddo enddo @@ -367,7 +367,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & else inci=ncells - incj=ncells + incj=ncells do j=1,nyc do i=1,nxc condition(i,j)=cape(inci/ncells,incj/ncells) @@ -390,7 +390,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif !nf - + !Vertical velocity has its own variable in order to condition on combination !of "condition" and vertical velocity. @@ -420,16 +420,16 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif !sgs/global -!Calculate neighbours and update the automata -!If ca-global is used, then nca independent CAs are called and weighted together to create one field; CA +!Calculate neighbours and update the automata +!If ca-global is used, then nca independent CAs are called and weighted together to create one field; CA call update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini,ilives, & nlives, ncells, nfracseed, nseed,nthresh, ca_global, & ca_sgs,nspinup, condition, vertvelhigh,nf,nca_plumes) - - if(ca_global == .true.)then + + if(ca_global)then CAstore(:,:) = CAstore(:,:) + CA(:,:) - elseif(ca_sgs == .true.)then + elseif(ca_sgs)then if(nf==1)then CA_DEEP(:,:)=CA(:,:) elseif(nf==2)then @@ -447,13 +447,13 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo !nf (nca) - if(ca_global == .true.)then + if(ca_global)then CAavg = CAstore / real(nca) endif !smooth CA field -if (ca_smooth ==.true. .and. ca_global ==.true.) then +if (ca_smooth .and. ca_global) then field_in=0. !get halo @@ -471,7 +471,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & do i=1,nlon ih=i+halo jh=j+halo - field_smooth(i,j)=(4.0*field_out(ih,jh,1)+2.0*field_out(ih-1,jh,1)+ & + field_smooth(i,j)=(4.0*field_out(ih,jh,1)+2.0*field_out(ih-1,jh,1)+ & 2.0*field_out(ih,jh-1,1)+2.0*field_out(ih+1,jh,1)+& 2.0*field_out(ih,jh+1,1)+2.0*field_out(ih-1,jh-1,1)+& 2.0*field_out(ih-1,jh+1,1)+2.0*field_out(ih+1,jh+1,1)+& @@ -527,7 +527,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !endif !Set the range for the nca individual ca_sgs patterns: -if(ca_sgs ==.true.)then +if(ca_sgs)then Detmax(1)=maxval(CA_DEEP(:,:)) call mp_reduce_max(Detmax(1)) @@ -535,7 +535,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon if(CA_DEEP(i,j)>0.)then - CA_DEEP(i,j)=CA_DEEP(i,j)/Detmax(1) !Now the range goes from 0-1 + CA_DEEP(i,j)=CA_DEEP(i,j)/Detmax(1) !Now the range goes from 0-1 endif enddo enddo @@ -560,7 +560,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon if(CA_DEEP(i,j)>0.)then - CA_DEEP(i,j)=(CA_DEEP(i,j)-CAmean) !Can we compute the median? + CA_DEEP(i,j)=(CA_DEEP(i,j)-CAmean) !Can we compute the median? endif enddo enddo @@ -568,8 +568,8 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !!! -!This is used for coupling with the Chikira-Sugiyama deep -!cumulus scheme. +!This is used for coupling with the Chikira-Sugiyama deep +!cumulus scheme. do j=1,nlat do i=1,nlon if(ca_plumes(i,j)==0)then @@ -582,7 +582,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & !Put back into blocks 1D array to be passed to physics !or diagnostics output - + do blk = 1, Atm_block%nblks do ix = 1,Atm_block%blksz(blk) i = Atm_block%index(blk)%ii(ix) - isc + 1 diff --git a/compns_stochy.F90 b/compns_stochy.F90 index af17eb0d..420ccdba 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -1,3 +1,9 @@ +module compns_stochy_mod + + implicit none + + contains + !----------------------------------------------------------------------- subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) !$$$ Subprogram Documentation Block @@ -27,12 +33,12 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! !$$$ - + use stochy_namelist_def - + implicit none - + integer, intent(out) :: iret integer, intent(in) :: nlunit,me,sz_nml character(len=*), intent(in) :: input_nml_file(sz_nml) @@ -44,7 +50,7 @@ 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,& + iseed_shum,iseed_sppt,shum,shum_tau,& shum_lscale,fhstoch,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,& @@ -62,7 +68,7 @@ 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 + ! mg, sfcperts pertz0 = -999. ! momentum roughness length amplitude pertshc = -999. ! soil hydraulic conductivity amp pertzt = -999. ! mom/heat roughness length amplitude @@ -199,3 +205,5 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! return end subroutine compns_stochy + +end module compns_stochy_mod diff --git a/dezouv_stochy.f b/dezouv_stochy.f index b95a3660..32f4e7db 100644 --- a/dezouv_stochy.f +++ b/dezouv_stochy.f @@ -1,10 +1,16 @@ + module dezouv_stochy_mod + + implicit none + + contains + subroutine dezouv_stochy(dev,zod,uev,vod,epsedn,epsodn, & snnp1ev,snnp1od,ls_node) cc cc use stochy_resol_def - use spectral_layout + use spectral_layout_mod use machine implicit none cc @@ -259,3 +265,5 @@ subroutine dezouv_stochy(dev,zod,uev,vod,epsedn,epsodn, cc return end + + end module dezouv_stochy_mod diff --git a/dozeuv_stochy.f b/dozeuv_stochy.f index deaf32a7..4ff5ad8f 100644 --- a/dozeuv_stochy.f +++ b/dozeuv_stochy.f @@ -1,8 +1,14 @@ + module dozeuv_stochy_mod + + implicit none + + contains + subroutine dozeuv_stochy(dod,zev,uod,vev,epsedn,epsodn, & snnp1ev,snnp1od,ls_node) cc use stochy_resol_def - use spectral_layout + use spectral_layout_mod use machine implicit none cc @@ -257,3 +263,5 @@ subroutine dozeuv_stochy(dod,zev,uod,vev,epsedn,epsodn, cc return end + + end module dozeuv_stochy_mod diff --git a/epslon_stochy.f b/epslon_stochy.f index 5ee6df2d..c7aace51 100644 --- a/epslon_stochy.f +++ b/epslon_stochy.f @@ -1,8 +1,14 @@ + module epslon_stochy_mod + + implicit none + + contains + subroutine epslon_stochy(epse,epso,epsedn,epsodn, & ls_node) cc use stochy_resol_def - use spectral_layout + use spectral_layout_mod use machine implicit none cc @@ -83,3 +89,5 @@ subroutine epslon_stochy(epse,epso,epsedn,epsodn, cc return end + + end module epslon_stochy_mod diff --git a/four_to_grid_stochy.f b/four_to_grid_stochy.F similarity index 75% rename from four_to_grid_stochy.f rename to four_to_grid_stochy.F index 8cf72290..ff3fdf38 100644 --- a/four_to_grid_stochy.f +++ b/four_to_grid_stochy.F @@ -1,4 +1,10 @@ + module four_to_grid_mod + use spectral_layout_mod, only: num_parthds_stochy => ompthreads + + implicit none + + contains subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, & lon_dim_coef,lon_dim_grid,lons_lat,lot) @@ -13,29 +19,45 @@ subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, integer lons_lat integer lot !________________________________________________________ +#ifdef MKL + integer*8 plan +#else real(kind=kind_dbl_prec) aux1crs(42002) real(kind=kind_dbl_prec) scale_ibm integer ibmsign integer init +#endif integer lot_thread integer num_threads integer nvar_thread_max integer nvar_1 integer nvar_2 integer thread - integer num_parthds_stochy +#ifdef MKL + include "fftw/fftw3.f" + integer NULL +#else + external dcrft + external scrft +#endif !________________________________________________________ - num_threads = min(num_parthds_stochy(),lot) - + num_threads = min(num_parthds_stochy,lot) + nvar_thread_max = (lot+num_threads-1)/num_threads - + if ( kind_dbl_prec == 8 ) then !------------------------------------ +#ifdef MKL +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,plan) +#else !$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) !$omp+shared(lon_dim_coef,lon_dim_grid) !$omp+shared(lot,num_threads,nvar_thread_max) !$omp+shared(ibmsign,scale_ibm) !$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) - +#endif do thread=1,num_threads ! start of thread loop .............. nvar_1=(thread-1)*nvar_thread_max + 1 nvar_2=min(nvar_1+nvar_thread_max-1,lot) @@ -43,6 +65,20 @@ subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, lot_thread=nvar_2 - nvar_1 + 1 if (nvar_2 >= nvar_1) then +#ifdef MKL + !call dfftw_plan_many_dft_c2r( + ! plan, 1, N, m, & + ! X, NULL, 1, dimx, & + ! Y, NULL, 1, dimy, & + ! fftw_flag) + call dfftw_plan_many_dft_c2r( & + & plan, 1, lons_lat, lot_thread, & + & syn_gr_a_1, NULL, 1, size(syn_gr_a_1,dim=1), & + & syn_gr_a_2, NULL, 1, size(syn_gr_a_2,dim=1), & + & FFTW_ESTIMATE) + call dfftw_execute(plan) + call dfftw_destroy_plan(plan) +#else init = 1 ibmsign = -1 scale_ibm = 1.0d0 @@ -61,23 +97,44 @@ subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, & lons_lat,lot_thread,ibmsign,scale_ibm, & aux1crs,22000, & aux1crs(22001),20000) +#endif endif enddo ! fin thread loop ...................................... else !------------------------------------------------------------ +#ifdef MKL +!$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) +!$omp+shared(lon_dim_coef,lon_dim_grid) +!$omp+shared(lot,num_threads,nvar_thread_max) +!$omp+private(thread,nvar_1,nvar_2,lot_thread,plan) +#else !$omp parallel do shared(syn_gr_a_1,syn_gr_a_2,lons_lat) !$omp+shared(lon_dim_coef,lon_dim_grid) !$omp+shared(lot,num_threads,nvar_thread_max) !$omp+shared(ibmsign,scale_ibm) !$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) - +#endif do thread=1,num_threads ! start of thread loop .............. nvar_1 = (thread-1)*nvar_thread_max + 1 nvar_2 = min(nvar_1+nvar_thread_max-1,lot) - if (nvar_2 >= nvar_1) then lot_thread = nvar_2 - nvar_1 + 1 - + + if (nvar_2 >= nvar_1) then +#ifdef MKL + !call sfftw_plan_many_dft_c2r( + ! plan, 1, N, m, & + ! X, NULL, 1, dimx, & + ! Y, NULL, 1, dimy, & + ! fftw_flag) + call sfftw_plan_many_dft_c2r( & + & plan, 1, lons_lat, lot_thread, & + & syn_gr_a_1, NULL, 1, size(syn_gr_a_1,dim=1), & + & syn_gr_a_2, NULL, 1, size(syn_gr_a_2,dim=1), & + & FFTW_ESTIMATE) + call sfftw_execute(plan) + call sfftw_destroy_plan(plan) +#else init = 1 ibmsign = -1 scale_ibm = 1.0d0 @@ -96,7 +153,7 @@ subroutine four_to_grid(syn_gr_a_1,syn_gr_a_2, & aux1crs,22000, & aux1crs(22001),20000, & aux1crs(22001),0) - +#endif endif enddo ! fin thread loop ...................................... endif !----------------------------------------------------------- @@ -125,26 +182,31 @@ subroutine grid_to_four(anl_gr_a_2,anl_gr_a_1, integer nvar_thread_max integer nvar_1,nvar_2 integer thread - integer num_parthds_stochy !________________________________________________________ - num_threads=min(num_parthds_stochy(),lot) - +#ifdef MKL + write(0,*) "ERROR in grid_to_four: srcft and drcft ", + & " must be replaced with MKL's FFTW calls. ABORT." + call sleep(5) + stop +#endif + num_threads=min(num_parthds_stochy,lot) + nvar_thread_max=(lot+num_threads-1)/num_threads - + if ( kind_dbl_prec == 8 ) then !------------------------------------ !$omp parallel do shared(anl_gr_a_1,anl_gr_a_2,lons_lat) !$omp+shared(lon_dim_coef,lon_dim_grid) !$omp+shared(lot,num_threads,nvar_thread_max) !$omp+shared(ibmsign,scale_ibm,rone) !$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) - + do thread=1,num_threads ! start of thread loop .............. nvar_1 = (thread-1)*nvar_thread_max + 1 nvar_2 = min(nvar_1+nvar_thread_max-1,lot) if (nvar_2 >= nvar_1) then lot_thread = nvar_2 - nvar_1 + 1 - + init = 1 ibmsign = 1 rone = 1.0d0 @@ -162,7 +224,7 @@ subroutine grid_to_four(anl_gr_a_2,anl_gr_a_1, & lons_lat,lot_thread,ibmsign,scale_ibm, & aux1crs,22000, & aux1crs(22001),20000) - + endif enddo ! fin thread loop ...................................... else !------------------------------------------------------------ @@ -171,14 +233,14 @@ subroutine grid_to_four(anl_gr_a_2,anl_gr_a_1, !$omp+shared(lot,num_threads,nvar_thread_max) !$omp+shared(ibmsign,scale_ibm,rone) !$omp+private(thread,nvar_1,nvar_2,lot_thread,init,aux1crs) - + do thread=1,num_threads ! start of thread loop .............. nvar_1 = (thread-1)*nvar_thread_max + 1 nvar_2 = min(nvar_1+nvar_thread_max-1,lot) if (nvar_2 >= nvar_1) then lot_thread=nvar_2 - nvar_1 + 1 - + init = 1 ibmsign = 1 rone = 1.0d0 @@ -198,10 +260,12 @@ subroutine grid_to_four(anl_gr_a_2,anl_gr_a_1, & aux1crs,22000, & aux1crs(22001),20000, & aux1crs(22001),0) - + endif enddo ! fin thread loop ...................................... endif !----------------------------------------------------------- !! return end + + end module four_to_grid_mod diff --git a/get_lats_node_a_stochy.f b/get_lats_node_a_stochy.f index b06c6419..86e85ec1 100644 --- a/get_lats_node_a_stochy.f +++ b/get_lats_node_a_stochy.f @@ -1,9 +1,15 @@ + module get_lats_node_a_stochy_mod + + implicit none + + contains + subroutine get_lats_node_a_stochy(me_fake,global_lats_a, & lats_nodes_a_fake,gl_lats_index, & global_time_sort_index,iprint) cc use stochy_resol_def - use spectral_layout + use spectral_layout_mod implicit none cc integer gl_lats_index,gl_start @@ -34,7 +40,7 @@ subroutine get_lats_node_a_stochy(me_fake,global_lats_a, lat = 1 nodes_tmp = nodes !jw if (liope .and. icolor .eq. 2) nodes_tmp = nodes -1 - + gl_start = gl_lats_index cc............................................. do ijk=1,latg @@ -72,7 +78,7 @@ subroutine get_lats_node_a_stochy(me_fake,global_lats_a, c$$$ . lats_nodes_a_fake endif enddo - + if(iprint.eq.1) print 220 220 format ('completed loop 200 in get_lats_a ') c @@ -82,3 +88,5 @@ subroutine get_lats_node_a_stochy(me_fake,global_lats_a, cc return end + + end module get_lats_node_a_stochy_mod diff --git a/get_ls_node_stochy.f b/get_ls_node_stochy.f index 82b396f6..51d9f85c 100644 --- a/get_ls_node_stochy.f +++ b/get_ls_node_stochy.f @@ -1,8 +1,14 @@ + module get_ls_node_stochy_mod + + implicit none + + contains + subroutine get_ls_node_stochy(me_fake,ls_node,ls_max_node_fake, c iprint) ! use stochy_resol_def - use spectral_layout + use spectral_layout_mod implicit none ! integer me_fake, ls_max_node_fake, iprint @@ -71,3 +77,5 @@ subroutine get_ls_node_stochy(me_fake,ls_node,ls_max_node_fake, ! return end + + end module get_ls_node_stochy_mod diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index f4f56237..d8c080fe 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -1,15 +1,25 @@ module get_stochy_pattern_mod - use stochy_resol_def - use spectral_layout - use stochy_namelist_def - use stochy_data_mod - use stochy_gg_def - use stochy_patterngenerator - use stochy_internal_state_mod + use machine, only : kind_dbl_prec, kind_evod + use stochy_resol_def, only : latg, latg2, levs, lonf, skeblevs + 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_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 + use stochy_gg_def, only : coslat_a + use stochy_patterngenerator_mod, only: random_pattern, ndimspec, & + patterngenerator_advance + use stochy_internal_state_mod, only: stochy_internal_state use fv_mp_mod, only : mp_reduce_sum,is_master - use fv_arrays_mod, only: fv_atmos_type use GFS_typedefs, only: GFS_control_type, GFS_grid_type use mersenne_twister, only: random_seed + use dezouv_stochy_mod, only: dezouv_stochy + use dozeuv_stochy_mod, only: dozeuv_stochy + use four_to_grid_mod, only: four_to_grid + use sumfln_stochy_mod, only: sumfln_stochy implicit none private @@ -17,7 +27,6 @@ module get_stochy_pattern_mod public get_random_pattern_sfc_fv3 public dump_patterns logical :: first_call=.true. -#include "mpif.h" contains subroutine get_random_pattern_fv3(rpattern,npatterns,& @@ -44,7 +53,6 @@ subroutine get_random_pattern_fv3(rpattern,npatterns,& real(kind=kind_dbl_prec) :: pattern_2d(nblks,maxlen) real(kind=kind_dbl_prec) :: pattern_1d(maxlen) real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: rslmsk - real(kind=kind_dbl_prec), allocatable, dimension(:) :: slmask,tlats,tlons integer :: blk kmsk0 = 0 @@ -67,7 +75,7 @@ subroutine get_random_pattern_fv3(rpattern,npatterns,& workg(i,lat) = glolal(i,j) enddo enddo - + call mp_reduce_sum(workg,lonf,latg) ! interpolate to cube grid @@ -76,19 +84,12 @@ subroutine get_random_pattern_fv3(rpattern,npatterns,& do blk=1,nblks len=size(Grid(blk)%xlat,1) pattern_1d = 0 - allocate(SLMASK(len)) - allocate(tlats(len)) - allocate(tlons(len)) - tlats=Grid(blk)%xlat*rad2deg - tlons=Grid(blk)%xlon*rad2deg - SLMASK = 0 - call la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) pattern_2d(blk,:)=pattern_1d(:) - deallocate(SLMASK) - deallocate(tlats) - deallocate(tlons) + end associate enddo deallocate(rslmsk) deallocate(workg) @@ -120,7 +121,6 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& 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 - real(kind=kind_dbl_prec), allocatable, dimension(:) :: slmask,tlats,tlons integer :: blk do k=1,nsfcpert @@ -154,19 +154,12 @@ subroutine get_random_pattern_sfc_fv3(rpattern,npatterns,& do blk=1,nblks len=size(Grid(blk)%xlat,1) pattern_1d = 0 - allocate(SLMASK(len)) - allocate(tlats(len)) - allocate(tlons(len)) - tlats=Grid(blk)%xlat*rad2deg - tlons=Grid(blk)%xlon*rad2deg - SLMASK = 0 - call la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) pattern_3d(blk,:,k)=pattern_1d(:) - deallocate(SLMASK) - deallocate(tlats) - deallocate(tlons) + 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)) deallocate(rslmsk) @@ -195,7 +188,6 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& 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 - real(kind=kind_dbl_prec), allocatable, dimension(:) :: slmask,tlats,tlons integer i,j,l,lat,ierr,n,nn,k,nt real(kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2du,wrk2dv @@ -245,23 +237,15 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& do blk=1,nblks len=size(Grid(blk)%xlat,1) pattern_1d = 0 - allocate(SLMASK(len)) - allocate(tlats(len)) - allocate(tlons(len)) - tlats=Grid(blk)%xlat*rad2deg - tlons=Grid(blk)%xlon*rad2deg - SLMASK = 0 - call la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) skebu_save(blk,:,k)=pattern_1d(:) - call la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) skebv_save(blk,:,k)=-1*pattern_1d(:) - deallocate(SLMASK) - deallocate(tlats) - deallocate(tlons) + end associate enddo enddo endif @@ -311,23 +295,15 @@ subroutine get_random_pattern_fv3_vect(rpattern,npatterns,& do blk=1,nblks len=size(Grid(blk)%xlat,1) pattern_1d = 0 - allocate(SLMASK(len)) - allocate(tlats(len)) - allocate(tlons(len)) - tlats=Grid(blk)%xlat*rad2deg - tlons=Grid(blk)%xlon*rad2deg - SLMASK = 0 - call la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + associate( tlats=>Grid(blk)%xlat*rad2deg,& + tlons=>Grid(blk)%xlon*rad2deg ) + call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) skebu_save(blk,:,skeblevs)=pattern_1d(:) - call la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,inttyp,& - pattern_1d(1:len),len,.false.,rslmsk,slmask,& - tlats,tlons,me) + call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,& + pattern_1d(1:len),len,rslmsk,tlats,tlons) skebv_save(blk,:,skeblevs)=-1*pattern_1d(:) - deallocate(SLMASK) - deallocate(tlats) - deallocate(tlons) + end associate enddo deallocate(rslmsk) deallocate(workgu) diff --git a/getcon_lag_stochy.f b/getcon_lag_stochy.f index fcab0b60..c92ed0d1 100644 --- a/getcon_lag_stochy.f +++ b/getcon_lag_stochy.f @@ -1,12 +1,20 @@ + module getcon_lag_stochy_mod + + implicit none + + contains + subroutine getcon_lag_stochy(lats_nodes_a,global_lats_a, & lats_nodes_h,global_lats_h_sn, & lonsperlat,xhalo,yhalo) use stochy_resol_def, only : jcap,latg,latg2,lonf - use spectral_layout, only : me,nodes + use spectral_layout_mod, only : me,nodes use stochy_gg_def, only : colrad_a,sinlat_a - use stochy_layout_lag, only : ipt_lats_node_h,lat1s_h,lats_dim_h, + use stochy_layout_lag, only : + & ipt_lats_node_h,lat1s_h,lats_dim_h, & lats_node_h,lats_node_h_max,lon_dim_h + use setlats_lag_stochy_mod, only: setlats_lag_stochy implicit none ! integer yhalo,xhalo @@ -77,3 +85,5 @@ subroutine getcon_lag_stochy(lats_nodes_a,global_lats_a, lon_dim_h = lonf + 1 + xhalo + xhalo !even/odd return end + + end module getcon_lag_stochy_mod diff --git a/getcon_spectral.F90 b/getcon_spectral.F90 index 85eb2223..387fc115 100644 --- a/getcon_spectral.F90 +++ b/getcon_spectral.F90 @@ -1,3 +1,9 @@ +module getcon_spectral_mod + + implicit none + + contains + subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & lats_nodes_a,global_lats_a, & lonsperlat,latsmax, & @@ -6,12 +12,19 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & snnp1ev,snnp1od, & plnev_a,plnod_a,pddev_a,pddod_a, & plnew_a,plnow_a,colat1) - + ! program log: ! 20110220 henry juang update code to fit mass_dp and ndslfv ! + use epslon_stochy_mod, only: epslon_stochy + use get_lats_node_a_stochy_mod, only: get_lats_node_a_stochy + use get_ls_node_stochy_mod, only: get_ls_node_stochy + use glats_stochy_mod, only: glats_stochy + use gozrineo_a_stochy_mod, only: gozrineo_a_stochy + use pln2eo_a_stochy_mod, only: pln2eo_a_stochy + use setlats_a_stochy_mod, only: setlats_a_stochy use stochy_resol_def - use spectral_layout + use spectral_layout_mod use stochy_gg_def use stochy_internal_state_mod @@ -32,7 +45,7 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & integer global_lats_ext(latg+2*jintmx+2*nypt*(nodes-1)) ! real(kind=kind_dbl_prec), dimension(len_trie_ls) :: epse, epsedn, snnp1ev - real(kind=kind_dbl_prec), dimension(len_trio_ls) :: epso, epsodn, snnp1od + real(kind=kind_dbl_prec), dimension(len_trio_ls) :: epso, epsodn, snnp1od ! real(kind=kind_dbl_prec), dimension(len_trie_ls,latg2) :: plnev_a, pddev_a, plnew_a real(kind=kind_dbl_prec), dimension(len_trio_ls,latg2) :: plnod_a, pddod_a, plnow_a @@ -69,7 +82,7 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & lonsperlat(latg+1-lat) = lonsperlat(lat) end do do node=1,nodes - call get_lats_node_a_stochy( node-1, global_lats_a,lats_nodes_a(node),& + call get_lats_node_a_stochy( node-1, global_lats_a,lats_nodes_a(node),& gl_lats_index,global_time_sort_index_a, iprint) enddo call setlats_a_stochy(lats_nodes_a,global_lats_a,iprint, lonsperlat) @@ -108,7 +121,7 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & lats_node_a_max = max(lats_node_a_max, lats_nodes_a(i)) enddo latsmax = lats_node_a_max - + ! ipt_lats_node_ext = 1 ! @@ -248,7 +261,7 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & 200 continue end do ! - + do j=1,lats_node_a lat = global_lats_a(ipt_lats_node_a-1+j) if ( lonsperlat(lat) == lonf ) then @@ -260,3 +273,5 @@ subroutine getcon_spectral ( ls_node,ls_nodes,max_ls_nodes, & ! return end + +end module getcon_spectral_mod diff --git a/glats_stochy.f b/glats_stochy.f index 113d0c44..4ed9d38f 100644 --- a/glats_stochy.f +++ b/glats_stochy.f @@ -1,3 +1,9 @@ + module glats_stochy_mod + + implicit none + + contains + subroutine glats_stochy(lgghaf,colrad,wgt,wgtcs,rcs2,iprint) ! ! Jan 2013 Henry Juang increase precision by kind_qdt_prec=16 @@ -100,3 +106,4 @@ subroutine poly(n,rad,p) return end + end module glats_stochy_mod diff --git a/gozrineo_stochy.f b/gozrineo_stochy.f index 79c4a212..01210ff9 100644 --- a/gozrineo_stochy.f +++ b/gozrineo_stochy.f @@ -1,10 +1,16 @@ + module gozrineo_a_stochy_mod + + implicit none + + contains + subroutine gozrineo_a_stochy(plnev_a,plnod_a, & pddev_a,pddod_a, & plnew_a,plnow_a, & epse,epso,rcs2_a,wgt_a,ls_node,num_lat) cc use stochy_resol_def - use spectral_layout + use spectral_layout_mod use machine implicit none cc @@ -170,3 +176,5 @@ subroutine gozrineo_a_stochy(plnev_a,plnod_a, cc return end + + end module gozrineo_a_stochy_mod diff --git a/initialize_spectral_mod.F90 b/initialize_spectral_mod.F90 index ca811eb6..78cff2f3 100644 --- a/initialize_spectral_mod.F90 +++ b/initialize_spectral_mod.F90 @@ -16,16 +16,19 @@ module initialize_spectral_mod !!uses: ! use machine - use spectral_layout, only : ipt_lats_node_a, lats_node_a_max,lon_dim_a,len_trie_ls,len_trio_ls & - ,nodes,ls_max_node,lats_dim_a,ls_dim,nodes_comp,lat1s_a + use spectral_layout_mod, only : ipt_lats_node_a, lats_node_a_max,lon_dim_a,len_trie_ls,len_trio_ls & + ,nodes,ls_max_node,lats_dim_a,ls_dim,lat1s_a use stochy_layout_lag, only : lat1s_h use stochy_internal_state_mod - use spectral_layout,only:lon_dims_a + use spectral_layout_mod,only:lon_dims_a, num_parthds_stochy => ompthreads use stochy_resol_def use stochy_namelist_def use fv_mp_mod, only : is_master use stochy_gg_def, only : wgt_a,sinlat_a,coslat_a,colrad_a,wgtcs_a,rcs2_a,lats_nodes_h,global_lats_h - use mpp_mod + use getcon_spectral_mod, only: getcon_spectral + use get_ls_node_stochy_mod, only: get_ls_node_stochy + use getcon_lag_stochy_mod, only: getcon_lag_stochy + !use mpp_mod #ifndef IBM USE omp_lib #endif @@ -58,26 +61,26 @@ subroutine initialize_spectral(gis_stochy, rc) ! print*,'before allocate lonsperlat,',& ! allocated(gis_stochy%lonsperlat),'latg=',latg ! - gis_stochy%nodes=mpp_npes() +! gis_stochy%nodes=mpp_npes() ! print*,'mpp_npes=',mpp_npes() nodes = gis_stochy%nodes npe_single_member = gis_stochy%npe_single_member lon_dim_a = lon_s + 2 jcap=ntrunc - jcap1 = jcap+1 - jcap2 = jcap+2 + jcap1 = jcap+1 + jcap2 = jcap+2 latg = lat_s - latg2 = latg/2 + latg2 = latg/2 lonf = lon_s - lnt = jcap2*jcap1/2 - lnuv = jcap2*jcap1 - lnt2 = lnt + lnt - lnt22 = lnt2 + 1 - lnte = (jcap2/2)*((jcap2/2)+1)-1 - lnto = (jcap2/2)*((jcap2/2)+1)-(jcap2/2) - lnted = lnte - lntod = lnto + lnt = jcap2*jcap1/2 + lnuv = jcap2*jcap1 + lnt2 = lnt + lnt + lnt22 = lnt2 + 1 + lnte = (jcap2/2)*((jcap2/2)+1)-1 + lnto = (jcap2/2)*((jcap2/2)+1)-(jcap2/2) + lnted = lnte + lntod = lnto gis_stochy%lnt2 = lnt2 @@ -88,12 +91,8 @@ subroutine initialize_spectral(gis_stochy, rc) allocate(wgtcs_a(latg2)) allocate(rcs2_a(latg2)) -!! create io communicator and comp communicator -!! - nodes_comp=nodes -! ! if (is_master()) then -! print*,'number of threads is',num_parthds_stochy() +! print*,'number of threads is',num_parthds_stochy ! print*,'number of mpi procs is',nodes ! endif ! @@ -101,7 +100,7 @@ subroutine initialize_spectral(gis_stochy, rc) ! print*,'allocating lonsperlat',latg allocate(gis_stochy%lonsperlat(latg)) ! print*,'size=',size(gis_stochy%lonsperlat) - + inquire (file="lonsperlat.dat", exist=file_exists) if ( .not. file_exists ) then @@ -110,12 +109,20 @@ subroutine initialize_spectral(gis_stochy, rc) else open (iunit,file='lonsperlat.dat',status='old',form='formatted', & action='read',iostat=iret) - if (iret /= 0) call mpp_error(FATAL,'error while reading lonsperlat.dat') + if (iret /= 0) then + write(0,*) 'error while reading lonsperlat.dat' + rc = 1 + return + end if rewind iunit read (iunit,*,iostat=iret) latghf,(gis_stochy%lonsperlat(i),i=1,latghf) if (latghf+latghf /= latg) then write(0,*)' latghf=',latghf,' not equal to latg/2=',latg/2 - if (iret /= 0) call mpp_error(FATAL,'lonsperlat file has wrong size') + if (iret /= 0) then + write(0,*) 'lonsperlat file has wrong size' + rc = 1 + return + end if endif do i=1,latghf gis_stochy%lonsperlat(latg-i+1) = gis_stochy%lonsperlat(i) @@ -212,7 +219,11 @@ subroutine initialize_spectral(gis_stochy, rc) allocate(coslat_a(latg)) allocate(lat1s_h(0:jcap)) ! - if(gis_stochy%iret/=0) call mpp_error(FATAL,'incompatible namelist - aborted in stochy') + if(gis_stochy%iret/=0) then + write(0,*) 'incompatible namelist - aborted in stochy' + rc = 1 + return + end if !! gis_stochy%lats_nodes_ext = 0 call getcon_spectral(gis_stochy%ls_node, gis_stochy%ls_nodes, & @@ -247,8 +258,9 @@ subroutine initialize_spectral(gis_stochy, rc) ! if (gis_stochy%me == 0) then ! print*, ' lats_dim_a=', lats_dim_a, ' lats_node_a=', gis_stochy%lats_node_a -! endif +! endif rc=0 end subroutine initialize_spectral + end module initialize_spectral_mod diff --git a/makefile b/makefile index 7d522de1..c8639292 100644 --- a/makefile +++ b/makefile @@ -4,11 +4,11 @@ inside_nems := $(wildcard ../../../conf/configure.nems) ifneq ($(strip $(inside_nems)),) include ../../../conf/configure.nems else - exist_configure_fv3 := $(wildcard ../conf/configure.fv3) + exist_configure_fv3 := $(wildcard ../FV3/conf/configure.fv3) ifneq ($(strip $(exist_configure_fv3)),) - include ../conf/configure.fv3 + include ../FV3/conf/configure.fv3 else - $(error "../conf/configure.fv3 file is missing. Run ./configure") + $(error "../FV3/conf/configure.fv3 file is missing. Run ./configure") endif $(info ) $(info Build standalone FV3 stochastic_physics ...) @@ -17,7 +17,7 @@ endif LIBRARY = libstochastic_physics.a -FFLAGS += -I./include -I./mpp/include -I../gfsphysics/ -I../atmos_cubed_sphere -I$(FMS_DIR) -I../namphysics +FFLAGS += -I./include -I../FV3/gfsphysics/ -I../FV3/atmos_cubed_sphere -I$(FMS_DIR) -I../FV3/namphysics SRCS_F = @@ -27,7 +27,6 @@ SRCS_f90 = \ ./plumes.f90 SRCS_f = \ - ./spectral_layout.f \ ./stochy_gg_def.f \ ./stochy_resol_def.f \ ./stochy_layout_lag.f \ @@ -44,9 +43,10 @@ SRCS_f = \ ./getcon_lag_stochy.f \ ./pln2eo_stochy.f \ ./dozeuv_stochy.f \ - ./dezouv_stochy.f + ./dezouv_stochy.f SRCS_F90 = \ + ./spectral_layout.F90 \ ./getcon_spectral.F90 \ ./stochy_namelist_def.F90 \ ./compns_stochy.F90 \ @@ -55,9 +55,9 @@ SRCS_F90 = \ ./stochy_patterngenerator.F90 \ ./stochy_data_mod.F90 \ ./get_stochy_pattern.F90 \ - ./initialize_spectral_mod.F90 + ./initialize_spectral_mod.F90 -SRCS_c = +SRCS_c = DEPEND_FILES = $(SRCS_f) $(SRCS_f90) $(SRCS_F) $(SRCS_F90) @@ -80,8 +80,8 @@ clean: @echo $(RM) -f $(LIBRARY) *__genmod.f90 *.o */*.o *.mod *.i90 *.lst *.i depend -MKDEPENDS = ../mkDepends.pl -include ../conf/make.rules +MKDEPENDS = ../FV3/mkDepends.pl +include ../FV3/conf/make.rules include ./depend diff --git a/num_parthds_stochy.f b/num_parthds_stochy.f index b6e6c295..2a71d453 100644 --- a/num_parthds_stochy.f +++ b/num_parthds_stochy.f @@ -6,5 +6,5 @@ function num_parthds_stochy() read(omp_threads,*,iostat=stat)number_of_openMP_threads num_parthds_stochy = number_of_openMP_threads return - end + end diff --git a/pln2eo_stochy.f b/pln2eo_stochy.f index bd254696..4c6576ba 100644 --- a/pln2eo_stochy.f +++ b/pln2eo_stochy.f @@ -1,3 +1,9 @@ + module pln2eo_a_stochy_mod + + implicit none + + contains + subroutine pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a, & ls_node,num_lat) ! @@ -5,7 +11,7 @@ subroutine pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a, ! underflow and overflow if necessary by henry juang 2012 july ! use stochy_resol_def - use spectral_layout + use spectral_layout_mod use machine implicit none ! @@ -277,3 +283,5 @@ subroutine pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a, cc return end + + end module pln2eo_a_stochy_mod diff --git a/plumes.f90 b/plumes.f90 index 697fec72..410cc47c 100644 --- a/plumes.f90 +++ b/plumes.f90 @@ -1,174 +1,174 @@ -subroutine plumes(V,L,a,row,col,kend) -implicit none - -!!January 2018 adapted from Mathworkds "ISLANDS" - -!! plumes finds all islands of four-connected elements in a matrix. -!! plumes returns a matrix the same size as the input matrix, with all of -!! the four-connected elements assigned an arbitrary island number. A second -!! return argument is an nx3 matrix. Each row of this matrix has -!! information about a particular island: the first column is the island -!! number which corresponds to the numbers in the first return argument, the -!! second column is the number of elements in the island, the third column -!! is the value of the elements in that island (the elements of the input -!! matrix). plumes will also return a binary matrix with ones in the -!! positions of the elements of the largest four-connected island. The -!! largest four-connected island is determined first by value; for example -!! if there are 2 islands each with 5 elements, one of which is made up of -!! 6's and the other of which is made up of 4's, the island with the 6's -!! will be returned. In the case where there are 2 islands with the same -!! number of elements and the same value, an arbitrary choice will be -!! returned. - -integer, intent(in) :: row,col,kend -integer, intent(in), dimension(row,col) :: a -integer, intent(out), dimension(kend) :: V,L -integer :: cnt,pp,mm,kk -integer :: i,j,cntr,gg,hh,ii,jj,idxx,IDX -integer, dimension(row,col) :: AG - - - -do j=1,col - do i=1,row - AG(i,j) = 0 !! This matrix will hold all of the islands found. - enddo -enddo - - -do j=1,kend -V(j)=0 -enddo - -L = V !! Hold the number of elements in each island. -cntr = 1 !! Label the individual islands by the order they are found. - -!! Notes on comments in the loops: CRNT is the element of the matrix we are -!! currently looking at. RGHT is the element to the immediate right of -!! CRNT. RUD is the element to the right and up from the CRNT. -!! The RUD element is directly above RGHT. - -do gg = 1,col-1 !! Look along the first row. - if (a(1,gg) == a(1,gg+1))then !! CRNT matches RGHT. - if (AG(1,gg) == 0)then !! CRNT does not have an island number. - AG(1,gg) = cntr !! Assign an island number to CRNT. - AG(1,gg+1) = cntr !! Assign an island number to RGHT. - V(cntr) = a(1,gg) !! Assign the value of the island. - L(cntr) = 2 !! Add to the island count. - cntr = cntr + 1 !! Increment the counter. - else !! CRNT does have an island number. - AG(1,gg+1) = AG(1,gg) !! Assign the island number to RGHT. - L(AG(1,gg)) = L(AG(1,gg)) + 1 !! Add to the island count. - endif - endif -enddo - - -do hh = 1,row-1 !! Look down the first column. - if (a(hh,1)==a(hh+1,1))then !! CRNT matches 'RGHT'. - if (AG(hh,1)==0)then !! CRNT does not have an island number. - AG(hh,1) = cntr !! Assign an island number to CRNT. - AG(hh+1,1) = cntr !! Assign an island number to 'RGHT'. - V(cntr) = a(hh,1) !! Assign the value of the island. - L(cntr) = 2 !! Add to the island count. - cntr = cntr + 1 - else !! CRNT does have an island number. - AG(hh+1,1) = AG(hh,1) !! Assign an island number to 'RGHT'. - L(AG(hh,1)) = L(AG(hh,1)) + 1 !! Add to the island count. - endif - endif -enddo - - -!! Now we can look at the rest of the matrix. -do ii = 2,row !! Start on the second row. - do jj = 1,col-1 !! Start on the first column. - if (a(ii,jj)==a(ii,jj+1))then !! CRNT matches RGHT. - if (a(ii,jj)==a(ii-1,jj+1))then !! CRNT matches RUD too. - if (AG(ii,jj)==0 .and. AG(ii-1,jj+1)==0)then !! Both aren't yet grouped. - AG(ii,jj) = cntr !! Give CRNT a new island number. - AG(ii-1,jj+1) = cntr !! Give RUD the new island num. - AG(ii,jj+1) = cntr !! Give RGHT the new island number. - V(cntr) = a(ii,jj) !! Store the value of the island. - L(cntr) = 3 !! Number of members in the new island. - cntr = cntr + 1 !! Increment the counter. - elseif (AG(ii-1,jj+1)==0)then !! RUD not yet been grouped, CRNT is. - AG(ii-1,jj+1) = AG(ii,jj) !! Give RUD CRNTs isl. num. - AG(ii,jj+1) = AG(ii,jj) !! And RGHT as well. - L(AG(ii,jj)) = L(AG(ii,jj)) + 2 !! Add to island size. - elseif (AG(ii,jj)==0)then !! CRNT not yet grouped, RUD is grouped. - AG(ii,jj) = AG(ii-1,jj+1) !! Give CRNT RUDs isl. num. - AG(ii,jj+1) = AG(ii-1,jj+1) !! And RGHT as well. - L(AG(ii-1,jj+1)) = L(AG(ii-1,jj+1)) + 2 !! Add to cnt. - else !! Both CRNT and RUD have been grouped: merge islands. - !! First decide which island has the least members. - if (L(AG(ii-1,jj+1)) L(idxx))then - goto 101 !! Stop search do old islnds. - endif -101 endif - enddo - if (cnt > L(idxx))then - goto 102 !! Stop search do old island numbers. -102 endif - enddo - endif - AG(ii,jj+1) = IDX !! Give RGHT the island number. - L(IDX) = L(IDX) + cnt !! Add to count. - L(idxx) = L(idxx) - cnt + 1 !! subtract from old island. - endif - elseif (AG(ii,jj)==0)then !! RGHT matches CRNT, not RUD. Need new isl. - AG(ii,jj) = cntr !! Give CRNT a new island number. - AG(ii,jj+1) = cntr !! Give RGHT the new island number. - V(cntr) = a(ii,jj) !! Store the new island number. - L(cntr) = 2 !! The number of members in the new island. - cntr = cntr + 1 !! Increment the counter. - else !! RGHT matches CRNT, not RUD. No new island needed. - AG(ii,jj+1) = AG(ii,jj) !! Give RGHT CRNTs island number. - L(AG(ii,jj)) = L(AG(ii,jj)) + 1 !! Add to island count. - endif - elseif (a(ii,jj+1)==a(ii-1,jj+1))then !! RUD & RGHT match, not CRNT. - if (AG(ii-1,jj+1)==0)then !! RUD has not yet been grouped. - AG(ii,jj+1) = cntr !! Give RGHT new island number. - AG(ii-1,jj+1) = cntr !! Give RUD new island number. - V(cntr) = a(ii,jj+1) !! Store the value of the island. - L(cntr) = 2 !! Add to island count. - cntr = cntr + 1 - else !! RUD is already part of a island. - AG(ii,jj+1) = AG(ii-1,jj+1) !! Add RGHT to RUD's island. - L(AG(ii-1,jj+1)) = L(AG(ii-1,jj+1)) + 1 !! Add island cnt. - endif - endif - enddo -enddo - - -! NSV = [(1:cntr-1)',L(1:cntr-1)',V(1:cntr-1)'] !!2nd output. -! NSV(NSV(:,2)==0,:)=[] !! Clear empty islands. - - -end subroutine plumes +subroutine plumes(V,L,a,row,col,kend) +implicit none + +!!January 2018 adapted from Mathworkds "ISLANDS" + +!! plumes finds all islands of four-connected elements in a matrix. +!! plumes returns a matrix the same size as the input matrix, with all of +!! the four-connected elements assigned an arbitrary island number. A second +!! return argument is an nx3 matrix. Each row of this matrix has +!! information about a particular island: the first column is the island +!! number which corresponds to the numbers in the first return argument, the +!! second column is the number of elements in the island, the third column +!! is the value of the elements in that island (the elements of the input +!! matrix). plumes will also return a binary matrix with ones in the +!! positions of the elements of the largest four-connected island. The +!! largest four-connected island is determined first by value; for example +!! if there are 2 islands each with 5 elements, one of which is made up of +!! 6's and the other of which is made up of 4's, the island with the 6's +!! will be returned. In the case where there are 2 islands with the same +!! number of elements and the same value, an arbitrary choice will be +!! returned. + +integer, intent(in) :: row,col,kend +integer, intent(in), dimension(row,col) :: a +integer, intent(out), dimension(kend) :: V,L +integer :: cnt,pp,mm,kk +integer :: i,j,cntr,gg,hh,ii,jj,idxx,IDX +integer, dimension(row,col) :: AG + + + +do j=1,col + do i=1,row + AG(i,j) = 0 !! This matrix will hold all of the islands found. + enddo +enddo + + +do j=1,kend +V(j)=0 +enddo + +L = V !! Hold the number of elements in each island. +cntr = 1 !! Label the individual islands by the order they are found. + +!! Notes on comments in the loops: CRNT is the element of the matrix we are +!! currently looking at. RGHT is the element to the immediate right of +!! CRNT. RUD is the element to the right and up from the CRNT. +!! The RUD element is directly above RGHT. + +do gg = 1,col-1 !! Look along the first row. + if (a(1,gg) == a(1,gg+1))then !! CRNT matches RGHT. + if (AG(1,gg) == 0)then !! CRNT does not have an island number. + AG(1,gg) = cntr !! Assign an island number to CRNT. + AG(1,gg+1) = cntr !! Assign an island number to RGHT. + V(cntr) = a(1,gg) !! Assign the value of the island. + L(cntr) = 2 !! Add to the island count. + cntr = cntr + 1 !! Increment the counter. + else !! CRNT does have an island number. + AG(1,gg+1) = AG(1,gg) !! Assign the island number to RGHT. + L(AG(1,gg)) = L(AG(1,gg)) + 1 !! Add to the island count. + endif + endif +enddo + + +do hh = 1,row-1 !! Look down the first column. + if (a(hh,1)==a(hh+1,1))then !! CRNT matches 'RGHT'. + if (AG(hh,1)==0)then !! CRNT does not have an island number. + AG(hh,1) = cntr !! Assign an island number to CRNT. + AG(hh+1,1) = cntr !! Assign an island number to 'RGHT'. + V(cntr) = a(hh,1) !! Assign the value of the island. + L(cntr) = 2 !! Add to the island count. + cntr = cntr + 1 + else !! CRNT does have an island number. + AG(hh+1,1) = AG(hh,1) !! Assign an island number to 'RGHT'. + L(AG(hh,1)) = L(AG(hh,1)) + 1 !! Add to the island count. + endif + endif +enddo + + +!! Now we can look at the rest of the matrix. +do ii = 2,row !! Start on the second row. + do jj = 1,col-1 !! Start on the first column. + if (a(ii,jj)==a(ii,jj+1))then !! CRNT matches RGHT. + if (a(ii,jj)==a(ii-1,jj+1))then !! CRNT matches RUD too. + if (AG(ii,jj)==0 .and. AG(ii-1,jj+1)==0)then !! Both aren't yet grouped. + AG(ii,jj) = cntr !! Give CRNT a new island number. + AG(ii-1,jj+1) = cntr !! Give RUD the new island num. + AG(ii,jj+1) = cntr !! Give RGHT the new island number. + V(cntr) = a(ii,jj) !! Store the value of the island. + L(cntr) = 3 !! Number of members in the new island. + cntr = cntr + 1 !! Increment the counter. + elseif (AG(ii-1,jj+1)==0)then !! RUD not yet been grouped, CRNT is. + AG(ii-1,jj+1) = AG(ii,jj) !! Give RUD CRNTs isl. num. + AG(ii,jj+1) = AG(ii,jj) !! And RGHT as well. + L(AG(ii,jj)) = L(AG(ii,jj)) + 2 !! Add to island size. + elseif (AG(ii,jj)==0)then !! CRNT not yet grouped, RUD is grouped. + AG(ii,jj) = AG(ii-1,jj+1) !! Give CRNT RUDs isl. num. + AG(ii,jj+1) = AG(ii-1,jj+1) !! And RGHT as well. + L(AG(ii-1,jj+1)) = L(AG(ii-1,jj+1)) + 2 !! Add to cnt. + else !! Both CRNT and RUD have been grouped: merge islands. + !! First decide which island has the least members. + if (L(AG(ii-1,jj+1)) L(idxx))then + goto 101 !! Stop search do old islnds. + endif +101 endif + enddo + if (cnt > L(idxx))then + goto 102 !! Stop search do old island numbers. +102 endif + enddo + endif + AG(ii,jj+1) = IDX !! Give RGHT the island number. + L(IDX) = L(IDX) + cnt !! Add to count. + L(idxx) = L(idxx) - cnt + 1 !! subtract from old island. + endif + elseif (AG(ii,jj)==0)then !! RGHT matches CRNT, not RUD. Need new isl. + AG(ii,jj) = cntr !! Give CRNT a new island number. + AG(ii,jj+1) = cntr !! Give RGHT the new island number. + V(cntr) = a(ii,jj) !! Store the new island number. + L(cntr) = 2 !! The number of members in the new island. + cntr = cntr + 1 !! Increment the counter. + else !! RGHT matches CRNT, not RUD. No new island needed. + AG(ii,jj+1) = AG(ii,jj) !! Give RGHT CRNTs island number. + L(AG(ii,jj)) = L(AG(ii,jj)) + 1 !! Add to island count. + endif + elseif (a(ii,jj+1)==a(ii-1,jj+1))then !! RUD & RGHT match, not CRNT. + if (AG(ii-1,jj+1)==0)then !! RUD has not yet been grouped. + AG(ii,jj+1) = cntr !! Give RGHT new island number. + AG(ii-1,jj+1) = cntr !! Give RUD new island number. + V(cntr) = a(ii,jj+1) !! Store the value of the island. + L(cntr) = 2 !! Add to island count. + cntr = cntr + 1 + else !! RUD is already part of a island. + AG(ii,jj+1) = AG(ii-1,jj+1) !! Add RGHT to RUD's island. + L(AG(ii-1,jj+1)) = L(AG(ii-1,jj+1)) + 1 !! Add island cnt. + endif + endif + enddo +enddo + + +! NSV = [(1:cntr-1)',L(1:cntr-1)',V(1:cntr-1)'] !!2nd output. +! NSV(NSV(:,2)==0,:)=[] !! Clear empty islands. + + +end subroutine plumes diff --git a/setlats_a_stochy.f b/setlats_a_stochy.f index e2d4085e..c58dc56c 100644 --- a/setlats_a_stochy.f +++ b/setlats_a_stochy.f @@ -1,8 +1,14 @@ + module setlats_a_stochy_mod + + implicit none + + contains + subroutine setlats_a_stochy(lats_nodes_a,global_lats_a, & iprint,lonsperlat) ! use stochy_resol_def , only : latg,lonf - use spectral_layout , only : nodes,me + use spectral_layout_mod , only : nodes,me ! implicit none ! @@ -130,7 +136,7 @@ subroutine setlats_a_stochy(lats_nodes_a,global_lats_a, lats_hold(jpt+1-i,ii) = latg+1-lats_hold(i,node) enddo enddo - + endif !! @@ -185,3 +191,5 @@ subroutine setlats_a_stochy(lats_nodes_a,global_lats_a, return end + + end module setlats_a_stochy_mod diff --git a/setlats_lag_stochy.f b/setlats_lag_stochy.f index fc926cbc..ff4e4f0a 100644 --- a/setlats_lag_stochy.f +++ b/setlats_lag_stochy.f @@ -1,8 +1,14 @@ + module setlats_lag_stochy_mod + + implicit none + + contains + subroutine setlats_lag_stochy(lats_nodes_a, global_lats_a, & lats_nodes_h, global_lats_h, yhalo) ! use stochy_resol_def, only : latg - use spectral_layout, only : me,nodes + use spectral_layout_mod, only : me,nodes implicit none ! integer yhalo @@ -117,3 +123,5 @@ subroutine setlats_lag_stochy(lats_nodes_a, global_lats_a, ! return end + + end module setlats_lag_stochy_mod diff --git a/spectral_layout.F90 b/spectral_layout.F90 new file mode 100644 index 00000000..247b69f3 --- /dev/null +++ b/spectral_layout.F90 @@ -0,0 +1,278 @@ +module spectral_layout_mod + + implicit none + +! +! program log: +! 20161011 philip pegion : make stochastic pattern generator standalone +! +! 20190503 dom heinzeller : add ompthreads and stochy_la2ga; todo: cleanup nodes, me, ... (defined multiple times in several files) +! + integer :: nodes, & + me, & + lon_dim_a, & + ls_dim, & + ls_max_node, & + lats_dim_a, & + lats_node_a, & + lats_node_a_max, & + ipt_lats_node_a, & + len_trie_ls, & + len_trio_ls, & + len_trie_ls_max, & + len_trio_ls_max, & + me_l_0, & + lats_dim_ext, & + lats_node_ext, & + ipt_lats_node_ext + + integer :: ompthreads = 1 +! + integer, allocatable :: lat1s_a(:), lon_dims_a(:),lon_dims_ext(:) + +contains + +! + ! interpolation from lat/lon or gaussian grid to other lat/lon grid + ! + subroutine stochy_la2ga(regin,imxin,jmxin,rinlon,rinlat,rlon,rlat, & + gauout,len,rslmsk, outlat, outlon) + use machine , only : kind_io8, kind_io4 + implicit none + ! interface variables + real (kind=kind_io8), intent(in) :: regin(imxin,jmxin) + integer, intent(in) :: imxin + integer, intent(in) :: jmxin + real (kind=kind_io8), intent(in) :: rinlon(imxin) + real (kind=kind_io8), intent(in) :: rinlat(jmxin) + real (kind=kind_io8), intent(in) :: rlon + real (kind=kind_io8), intent(in) :: rlat + real (kind=kind_io8), intent(out) :: gauout(len) + integer, intent(in) :: len + real (kind=kind_io8), intent(in) :: rslmsk(imxin,jmxin) + real (kind=kind_io8), intent(in) :: outlat(len) + real (kind=kind_io8), intent(in) :: outlon(len) + ! local variables + real (kind=kind_io8) :: wei4,wei3,wei2,sum2,sum1,sum3,wei1,sum4 + real (kind=kind_io8) :: wsum,wsumiv,sums,sumn,wi2j2,x,y,wi1j1 + real (kind=kind_io8) :: wi1j2,wi2j1,aphi,rnume,alamd,denom + integer :: jy,ix,i,j,jq,jx + integer :: j1,j2,ii,i1,i2,kmami,it + integer :: nx,kxs,kxt + integer :: iindx1(len) + integer :: iindx2(len) + integer :: jindx1(len) + integer :: jindx2(len) + real(kind=kind_io8) :: ddx(len) + real(kind=kind_io8) :: ddy(len) + real(kind=kind_io8) :: wrk(len) + integer :: len_thread_m + integer :: len_thread + integer :: i1_t + integer :: i2_t +! + len_thread_m = (len+ompthreads-1) / ompthreads +! + !$omp parallel do num_threads(ompthreads) default(none) & + !$omp private(i1_t,i2_t,len_thread,it,i,ii,i1,i2) & + !$omp private(j,j1,j2,jq,ix,jy,nx,kxs,kxt,kmami) & + !$omp private(alamd,denom,rnume,aphi,x,y,wsum,wsumiv,sum1,sum2) & + !$omp private(sum3,sum4,wi1j1,wi2j1,wi1j2,wi2j2,wei1,wei2,wei3,wei4) & + !$omp private(sumn,sums) & + !$omp shared(imxin,jmxin) & + !$omp shared(outlon,outlat,wrk,iindx1,rinlon,jindx1,rinlat,ddx,ddy) & + !$omp shared(rlon,rlat,regin,gauout) & + !$omp shared(ompthreads,len_thread_m,len,iindx2,jindx2,rslmsk) + do it=1,ompthreads ! start of threaded loop + i1_t = (it-1)*len_thread_m+1 + i2_t = min(i1_t+len_thread_m-1,len) + len_thread = i2_t-i1_t+1 +! +! find i-index for interpolation +! + do i=i1_t, i2_t + alamd = outlon(i) + if (alamd .lt. rlon) alamd = alamd + 360.0 + if (alamd .gt. 360.0+rlon) alamd = alamd - 360.0 + wrk(i) = alamd + iindx1(i) = imxin + enddo + do i=i1_t,i2_t + do ii=1,imxin + if(wrk(i) .ge. rinlon(ii)) iindx1(i) = ii + enddo + enddo + do i=i1_t,i2_t + i1 = iindx1(i) + if (i1 .lt. 1) i1 = imxin + i2 = i1 + 1 + if (i2 .gt. imxin) i2 = 1 + iindx1(i) = i1 + iindx2(i) = i2 + denom = rinlon(i2) - rinlon(i1) + if(denom.lt.0.) denom = denom + 360. + rnume = wrk(i) - rinlon(i1) + if(rnume.lt.0.) rnume = rnume + 360. + ddx(i) = rnume / denom + enddo +! +! find j-index for interplation +! + if(rlat.gt.0.) then + do j=i1_t,i2_t + jindx1(j)=0 + enddo + do jx=1,jmxin + do j=i1_t,i2_t + if(outlat(j).le.rinlat(jx)) jindx1(j) = jx + enddo + enddo + do j=i1_t,i2_t + jq = jindx1(j) + aphi=outlat(j) + if(jq.ge.1 .and. jq .lt. jmxin) then + j2=jq+1 + j1=jq + ddy(j)=(aphi-rinlat(j1))/(rinlat(j2)-rinlat(j1)) + elseif (jq .eq. 0) then + j2=1 + j1=1 + if(abs(90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + else + j2=jmxin + j1=jmxin + if(abs(-90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(-90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + endif + jindx1(j)=j1 + jindx2(j)=j2 + enddo + else + do j=i1_t,i2_t + jindx1(j) = jmxin+1 + enddo + do jx=jmxin,1,-1 + do j=i1_t,i2_t + if(outlat(j).le.rinlat(jx)) jindx1(j) = jx + enddo + enddo + do j=i1_t,i2_t + jq = jindx1(j) + aphi=outlat(j) + if(jq.gt.1 .and. jq .le. jmxin) then + j2=jq + j1=jq-1 + ddy(j)=(aphi-rinlat(j1))/(rinlat(j2)-rinlat(j1)) + elseif (jq .eq. 1) then + j2=1 + j1=1 + if(abs(-90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(-90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + else + j2=jmxin + j1=jmxin + if(abs(90.-rinlat(j1)).gt.0.001) then + ddy(j)=(aphi-rinlat(j1))/(90.-rinlat(j1)) + else + ddy(j)=0.0 + endif + endif + jindx1(j)=j1 + jindx2(j)=j2 + enddo + endif +! + sum1 = 0. + sum2 = 0. + sum3 = 0. + sum4 = 0. + do i=1,imxin + sum1 = sum1 + regin(i,1) + sum2 = sum2 + regin(i,jmxin) + enddo + sum1 = sum1 / imxin + sum2 = sum2 / imxin + sum3 = sum1 + sum4 = sum2 +! +! quasi-bilinear interpolation +! + do i=i1_t,i2_t + y = ddy(i) + j1 = jindx1(i) + j2 = jindx2(i) + x = ddx(i) + i1 = iindx1(i) + i2 = iindx2(i) +! + wi1j1 = (1.-x) * (1.-y) + wi2j1 = x *( 1.-y) + wi1j2 = (1.-x) * y + wi2j2 = x * y +! + wsum = wi1j1 + wi2j1 + wi1j2 + wi2j2 + wrk(i) = wsum + if(wsum.ne.0.) then + wsumiv = 1./wsum + if(j1.ne.j2) then + gauout(i) = (wi1j1*regin(i1,j1) + wi2j1*regin(i2,j1) + & + wi1j2*regin(i1,j2) + wi2j2*regin(i2,j2)) & + *wsumiv + else + if (rlat .gt. 0.0) then + sumn = sum3 + sums = sum4 + if( j1 .eq. 1) then + gauout(i) = (wi1j1*sumn +wi2j1*sumn + & + wi1j2*regin(i1,j2)+wi2j2*regin(i2,j2)) & + * wsumiv + elseif (j1 .eq. jmxin) then + gauout(i) = (wi1j1*regin(i1,j1)+wi2j1*regin(i2,j1)+ & + wi1j2*sums +wi2j2*sums ) & + * wsumiv + endif + else + sums = sum3 + sumn = sum4 + if( j1 .eq. 1) then + gauout(i) = (wi1j1*regin(i1,j1)+wi2j1*regin(i2,j1)+ & + wi1j2*sums +wi2j2*sums ) & + * wsumiv + elseif (j1 .eq. jmxin) then + gauout(i) = (wi1j1*sumn +wi2j1*sumn + & + wi1j2*regin(i1,j2)+wi2j2*regin(i2,j2)) & + * wsumiv + endif + endif + endif ! if j1 .ne. j2 + endif + enddo + do i=i1_t,i2_t + j1 = jindx1(i) + j2 = jindx2(i) + i1 = iindx1(i) + i2 = iindx2(i) + if(wrk(i) .eq. 0.0) then + write(6,*) ' la2ga: bad rslmsk given' + call sleep(2) + stop + endif + enddo + enddo ! end of threaded loop +!$omp end parallel do +! + return +! + end subroutine stochy_la2ga + +end module spectral_layout_mod diff --git a/spectral_layout.f b/spectral_layout.f deleted file mode 100644 index bf064db0..00000000 --- a/spectral_layout.f +++ /dev/null @@ -1,28 +0,0 @@ - module spectral_layout - implicit none -! -! program log: -! 20161011 philip pegion : make stochastic pattern generator standalone -! - - integer nodes, nodes_comp,nodes_io, - & me,lon_dim_a, - & ls_dim, - & ls_max_node, - & lats_dim_a, - & lats_node_a, - & lats_node_a_max, - & ipt_lats_node_a, - & len_trie_ls, - & len_trio_ls, - & len_trie_ls_max, - & len_trio_ls_max, - & me_l_0, - - & lats_dim_ext, - & lats_node_ext, - & ipt_lats_node_ext -! - INTEGER ,ALLOCATABLE :: lat1s_a(:), lon_dims_a(:),lon_dims_ext(:) - - end module spectral_layout diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index a025a19a..bd0b627d 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -1,48 +1,79 @@ +module stochastic_physics -subroutine init_stochastic_physics(Model,Init_parm,nblks,Grid) +implicit none + +private + +public :: init_stochastic_physics, run_stochastic_physics + +contains + +subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) use fv_mp_mod, only : is_master 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,& rad2deg,INTTYP,wlon,rnlat,gis_stochy,vfact_skeb,vfact_sppt,vfact_shum,skeb_vpts,skeb_vwts,sl - -use get_stochy_pattern_mod,only : get_random_pattern_fv3_vect, get_random_pattern_sfc_fv3 use stochy_resol_def , only : latg,lonf,skeblevs use stochy_gg_def,only : colrad_a use stochy_namelist_def use physcons, only: con_pi -use spectral_layout,only:me +use spectral_layout_mod,only:me,ompthreads use mpp_mod -use MPI -use GFS_typedefs, only: GFS_control_type, GFS_init_type, GFS_coupling_type, GFS_grid_type +use GFS_typedefs, only: GFS_control_type, GFS_init_type implicit none type(GFS_control_type), intent(inout) :: Model -type(GFS_init_type), intent(in) :: Init_parm -integer,intent(in) :: nblks -type(GFS_grid_type), intent(in) :: Grid(nblks) +type(GFS_init_type), intent(in) :: Init_parm +integer, intent(in) :: ntasks +integer, intent(in) :: nthreads + +integer :: nblks +integer :: iret real*8 :: PRSI(Model%levs),PRSL(Model%levs),dx real, allocatable :: skeb_vloc(:) integer :: k,kflip,latghf,nodes,blk,k2 character*2::proc +! Set/update shared variables in spectral_layout_mod +ompthreads = nthreads + +! ------------------------------------------ + +nblks = size(Model%blksz) + ! replace rad2deg=180.0/con_pi INTTYP=0 ! bilinear interpolation -me=mpp_pe() -nodes=mpp_npes() +me=Model%me +nodes=ntasks gis_stochy%me=me gis_stochy%nodes=nodes -call init_stochdata(Model%levs,Model%dtp,Model%input_nml_file,Model%fn_nml,Init_parm%nlunit) +call init_stochdata(Model%levs,Model%dtp,Model%input_nml_file,Model%fn_nml,Init_parm%nlunit,iret) ! check to see decomposition !if(Model%isppt_deep == .true.)then !do_sppt = .true. !endif -Model%do_sppt=do_sppt +! check namelist entries for consistency +if (Model%do_sppt.neqv.do_sppt) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_sppt and sppt' + return +else if (Model%do_shum.neqv.do_shum) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_shum and shum' + return +else if (Model%do_skeb.neqv.do_skeb) then + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_skeb and skeb' + return +else if (Model%do_sfcperts.neqv.do_sfcperts) then ! mg, sfc-perts + write(0,'(*(a))') 'Logic error in stochastic_physics_init: incompatible', & + & ' namelist settings do_sfcperts and pertz0 / pertshc / pertzt / pertlai / pertvegf / pertalb' + return +end if +! update remaining model configuration parameters from namelist Model%use_zmtnblck=use_zmtnblck -Model%do_shum=do_shum -Model%do_skeb=do_skeb Model%skeb_npass=skeb_npass -Model%do_sfcperts=do_sfcperts ! mg, sfc-perts Model%nsfcpert=nsfcpert ! mg, sfc-perts Model%pertz0=pertz0 ! mg, sfc-perts Model%pertzt=pertzt ! mg, sfc-perts @@ -91,7 +122,7 @@ subroutine init_stochastic_physics(Model,Init_parm,nblks,Grid) else vfact_skeb(k) = 1.0 endif - if (is_master()) print *,'skeb vert profile',k,sl(k),vfact_skeb(k) + if (is_master()) print *,'skeb vert profile',k,sl(k),vfact_skeb(k) enddo ! calculate vertical interpolation weights do k=1,skeblevs @@ -110,8 +141,8 @@ subroutine init_stochastic_physics(Model,Init_parm,nblks,Grid) skeb_vpts(k,1)=k2 skeb_vwts(k,2)=(skeb_vloc(k2)-sl(k))/(skeb_vloc(k2)-skeb_vloc(k2+1)) ENDIF - ENDDO -ENDDO + ENDDO +ENDDO deallocate(skeb_vloc) if (is_master()) then DO k=1,Model%levs @@ -157,39 +188,38 @@ subroutine init_stochastic_physics(Model,Init_parm,nblks,Grid) end subroutine init_stochastic_physics -subroutine run_stochastic_physics(nblks,Model,Grid,Coupling) +subroutine run_stochastic_physics(Model, Grid, Coupling, nthreads) use fv_mp_mod, only : is_master 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,& - get_random_pattern_sfc_fv3 ! mg, sfc-perts +use get_stochy_pattern_mod,only : get_random_pattern_fv3,get_random_pattern_fv3_vect,dump_patterns use stochy_resol_def , only : latg,lonf use stochy_namelist_def -use spectral_layout,only:me +use spectral_layout_mod,only:me,ompthreads use mpp_mod -use MPI -use GFS_typedefs, only: GFS_control_type, GFS_grid_type,GFS_Coupling_type +use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_Coupling_type implicit none -integer,intent(in) :: nblks type(GFS_control_type), intent(in) :: Model -type(GFS_grid_type), intent(in) :: Grid(nblks) -type(GFS_coupling_type), intent(inout) :: Coupling(nblks) +type(GFS_grid_type), intent(in) :: Grid(:) +type(GFS_coupling_type), intent(inout) :: Coupling(:) +integer, intent(in) :: nthreads real,allocatable :: tmp_wts(:,:),tmpu_wts(:,:,:),tmpv_wts(:,:,:) !D-grid integer :: k integer j,ierr,i -integer :: blk, len, maxlen +integer :: nblks, blk, len, maxlen character*120 :: sfile character*6 :: STRFH if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return -maxlen = 0 -DO blk=1,nblks - maxlen = max(maxlen,size(Grid(blk)%xlat,1)) -ENDDO + +! Update number of threads in shared variables in spectral_layout_mod and set block-related variables +ompthreads = nthreads +nblks = size(Model%blksz) +maxlen = maxval(Model%blksz(:)) + ! check to see if it is time to write out random patterns if (Model%phour .EQ. fhstoch) then write(STRFH,FMT='(I6.6)') nint(Model%phour) @@ -235,35 +265,46 @@ subroutine run_stochastic_physics(nblks,Model,Grid,Coupling) end subroutine run_stochastic_physics -subroutine run_stochastic_physics_sfc(nblks,Model,Grid,Coupling) +end module stochastic_physics + + +module stochastic_physics_sfc + +implicit none + +private + +public :: run_stochastic_physics_sfc + +contains + +subroutine run_stochastic_physics_sfc(Model, Grid, Coupling) use fv_mp_mod, 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 -use spectral_layout,only:me -use mpp_mod -use MPI -use GFS_typedefs, only: GFS_control_type, GFS_grid_type,GFS_Coupling_type +!use mpp_mod +use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_Coupling_type implicit none -integer,intent(in) :: nblks type(GFS_control_type), intent(in) :: Model -type(GFS_grid_type), intent(in) :: Grid(nblks) -type(GFS_coupling_type), intent(inout) :: Coupling(nblks) +type(GFS_grid_type), intent(in) :: Grid(:) +type(GFS_coupling_type), intent(inout) :: Coupling(:) real,allocatable :: tmpsfc_wts(:,:,:) !D-grid integer :: k integer j,ierr,i -integer :: blk, len, maxlen +integer :: nblks, blk, len, maxlen character*120 :: sfile character*6 :: STRFH if (.NOT. do_sfcperts) return -maxlen = 0 -DO blk=1,nblks - maxlen = max(maxlen,size(Grid(blk)%xlat,1)) -ENDDO + +! Set block-related variables +nblks = size(Model%blksz) +maxlen = maxval(Model%blksz(:)) + allocate(tmpsfc_wts(nblks,maxlen,Model%nsfcpert)) ! mg, sfc-perts if (is_master()) then print*,'In init_stochastic_physics: do_sfcperts ',do_sfcperts @@ -282,3 +323,5 @@ subroutine run_stochastic_physics_sfc(nblks,Model,Grid,Coupling) 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 dca8762f..f1cdf87b 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -2,17 +2,19 @@ module stochy_data_mod ! set up and initialize stochastic random patterns. - use spectral_layout + use spectral_layout_mod, only: len_trie_ls,len_trio_ls,ls_dim,ls_max_node use stochy_resol_def, only : skeblevs,levs,jcap,lonf,latg use stochy_namelist_def use constants_mod, only : radius - use fv_mp_mod,only:mp_bcst - use stochy_patterngenerator, only: random_pattern, patterngenerator_init,& + use spectral_layout_mod, only : me, nodes + use fv_mp_mod, only: mp_bcst, is_master + use stochy_patterngenerator_mod, only: random_pattern, patterngenerator_init,& getnoise, patterngenerator_advance,ndimspec,chgres_pattern,computevarspec_r - use initialize_spectral_mod + use initialize_spectral_mod, only: initialize_spectral use stochy_internal_state_mod ! use mersenne_twister_stochy, only : random_seed use mersenne_twister, only : random_seed + use compns_stochy_mod, only : compns_stochy implicit none private @@ -33,9 +35,9 @@ module stochy_data_mod real(kind=kind_dbl_prec),public, allocatable :: skebu_save(:,:,:),skebv_save(:,:,:) integer,public :: INTTYP type(stochy_internal_state),public :: gis_stochy - + contains - subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) + subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) ! initialize random patterns. A spinup period of spinup_efolds times the ! temporal time scale is run for each pattern. @@ -43,12 +45,13 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) character(len=*), intent(in) :: input_nml_file(:) character(len=64), intent(in) :: fn_nml real, intent(in) :: delt + integer, intent(out) :: iret real :: pertsfc(1) real :: rnn1 - integer :: nn,nspinup,k,nm,spinup_efolds,stochlun,ierr,iret,n + integer :: nn,nspinup,k,nm,spinup_efolds,stochlun,ierr,n integer :: locl,indev,indod,indlsod,indlsev - integer :: l,jbasev,jbasod,nodes + integer :: l,jbasev,jbasod real(kind_dbl_prec),allocatable :: noise_e(:,:),noise_o(:,:) include 'function_indlsod' include 'function_indlsev' @@ -59,7 +62,6 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) if(is_master()) print*,'in init stochdata' call compns_stochy (me,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return - nodes=mpp_npes() if (nodes.GE.lat_s/2) then lat_s=(int(nodes/12)+1)*24 lon_s=lat_s*2 @@ -67,6 +69,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) if (is_master()) print*,'WARNING: spectral resolution is too low for number of mpi_tasks, resetting lon_s,lat_s,and ntrunc to',lon_s,lat_s,ntrunc endif call initialize_spectral(gis_stochy, iret) + if (iret/=0) return allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) ! determine number of random patterns to be used for each scheme. do n=1,size(sppt) @@ -125,7 +128,11 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) if (stochini) then print*,'opening stoch_ini' OPEN(stochlun,file='stoch_ini',form='unformatted',iostat=ierr,status='old') - if (ierr .NE. 0) call mpp_error(FATAL,'error opening stoch_ini') + if (ierr .NE. 0) then + write(0,*) 'error opening stoch_ini' + iret = ierr + return + end if endif endif ! no spinup needed if initial patterns are defined correctly. @@ -318,7 +325,7 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit) 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)) + 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 diff --git a/stochy_gg_def.f b/stochy_gg_def.f index 9470c436..912874e8 100644 --- a/stochy_gg_def.f +++ b/stochy_gg_def.f @@ -1,7 +1,7 @@ module stochy_gg_def use machine implicit none - + real(kind=kind_dbl_prec), allocatable, dimension(:) :: colrad_a, & wgt_a, wgtcs_a, rcs2_a, sinlat_a, coslat_a ! diff --git a/stochy_internal_state_mod.F90 b/stochy_internal_state_mod.F90 index 85a98da0..1b2f1cbb 100644 --- a/stochy_internal_state_mod.F90 +++ b/stochy_internal_state_mod.F90 @@ -1,6 +1,6 @@ ! -! !module: stochy_internal_state_mod +! !module: stochy_internal_state_mod ! --- internal state definition of the ! gridded component of the spectral random patterns ! @@ -13,15 +13,15 @@ ! ! !interface: ! - + module stochy_internal_state_mod !!uses: !------ - use spectral_layout + use spectral_layout_mod use stochy_gg_def use stochy_resol_def - + implicit none private @@ -120,7 +120,7 @@ module stochy_internal_state_mod integer ikey,nrank_all,kcolor real(kind=kind_dbl_prec) cons0p5,cons1200,cons3600,cons0 - + ! ! ----------------------------------------------------- end type stochy_internal_state ! end type define diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index 06fac4f4..230145fb 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -5,7 +5,7 @@ module stochy_namelist_def ! use machine implicit none - + public integer nsskeb,lon_s,lat_s,ntrunc diff --git a/stochy_patterngenerator.F90 b/stochy_patterngenerator.F90 index 2ade7b90..6c400af2 100644 --- a/stochy_patterngenerator.F90 +++ b/stochy_patterngenerator.F90 @@ -1,12 +1,12 @@ -module stochy_patterngenerator +module stochy_patterngenerator_mod ! generate random patterns with specified temporal and spatial auto-correlation ! in spherical harmonic space. use machine - use spectral_layout, only: len_trie_ls, len_trio_ls, ls_dim, ls_max_node + use spectral_layout_mod, only: len_trie_ls, len_trio_ls, ls_dim, ls_max_node ! use mersenne_twister_stochy, only: random_setseed,random_gauss,random_stat use mersenne_twister, only: random_setseed,random_gauss,random_stat - use fv_mp_mod + use fv_mp_mod,only: is_master, mp_bcst implicit none private @@ -31,7 +31,7 @@ module stochy_patterngenerator end type random_pattern integer :: nlons,nlats,ntrunc,ndimspec - + contains subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& @@ -54,16 +54,16 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& nlons = nlon nlats = nlat ntrunc = jcap - ndimspec = (ntrunc+1)*(ntrunc+2)/2 + ndimspec = (ntrunc+1)*(ntrunc+2)/2 ! propagate seed supplied from namelist to all patterns... if (iseed(1) .NE. 0) then do np=2,npatterns if (iseed(np).EQ.0) then - iseed(np)=iseed(1)+np*100000000 + iseed(np)=iseed(1)+np*100000000 endif enddo endif - + do np=1,npatterns allocate(rpattern(np)%idx(0:ntrunc,0:ntrunc)) allocate(rpattern(np)%idx_e(len_trie_ls)) @@ -247,8 +247,13 @@ subroutine getnoise(rpattern,noise_e,noise_o) endif enddo end subroutine getnoise - + subroutine patterngenerator_advance(rpattern,k,skeb_first_call) + +#ifdef TRANSITION +!DIR$ OPTIMIZE:1 +#endif + ! advance 1st-order autoregressive process with ! specified autocorrelation (phi) and variance spectrum (spectrum) real(kind_dbl_prec) :: noise_e(len_trie_ls,2) @@ -354,4 +359,4 @@ subroutine chgres_pattern(pattern2din,pattern2dout,ntruncin,ntruncout) deallocate(idxin) end subroutine chgres_pattern -end module stochy_patterngenerator +end module stochy_patterngenerator_mod diff --git a/stochy_resol_def.f b/stochy_resol_def.f index 708e31c8..efbbd42e 100644 --- a/stochy_resol_def.f +++ b/stochy_resol_def.f @@ -37,7 +37,7 @@ module stochy_resol_def ! integer kdpphi, kzzphi, kdplam, kzzlam ! - integer kdtphi, kdrphi, kdtlam, kdrlam + integer kdtphi, kdrphi, kdtlam, kdrlam integer kdulam, kdvlam, kduphi, kdvphi diff --git a/sumfln_stochy.f b/sumfln_stochy.f index ad6b03b0..91a751fa 100644 --- a/sumfln_stochy.f +++ b/sumfln_stochy.f @@ -1,3 +1,9 @@ + module sumfln_stochy_mod + + implicit none + + contains + subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, & nvars,ls_node,latl2, & workdim,nvarsdim,four_gr, @@ -7,16 +13,17 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, & lons_lat,londi,latl,nvars_0) ! use stochy_resol_def , only : jcap,latgd - use spectral_layout , only : len_trie_ls,len_trio_ls, + use spectral_layout_mod , only : len_trie_ls,len_trio_ls, & ls_dim,ls_max_node,me,nodes use machine - use fv_mp_mod - use mpp_mod - use fv_arrays_mod, only: fv_atmos_type + use spectral_layout_mod, only : num_parthds_stochy => ompthreads + !or : use fv_mp_mod ? + use mpp_mod, only: mpp_npes, mpp_alltoall, mpp_get_current_pelist implicit none ! - include 'mpif.h' + external esmf_dgemm +! integer lat1s(0:jcap),latl2 ! integer nvars,nvars_0 @@ -79,12 +86,11 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, include 'function_indlsod' ! real(kind=kind_dbl_prec), parameter :: cons0=0.0d0, cons1=1.0d0 - integer num_parthds_stochy ! ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! arrsz=2*nvars*ls_dim*workdim*nodes - num_threads = min(NUM_PARTHDS_STOCHY(),nvars) + num_threads = min(num_parthds_stochy,nvars) nvar_thread_max = (nvars+num_threads-1)/num_threads npes = mpp_npes() allocate(pelist(0:npes-1)) @@ -181,7 +187,7 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, do node = 1, nodes - 1 ilat_list(node+1) = ilat_list(node) + lats_nodes(node) end do - + !$omp parallel do private(node,jj,ilat,lat,ipt_ls,nvar,kn,n2) do node=1,nodes do jj=1,lats_nodes(node) @@ -240,7 +246,7 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, end do work1dr(1:arrsz)=>workr work1ds(1:arrsz)=>works - call MPP_ALLTOALL(work1ds, sendcounts, sdispls, + call mpp_alltoall(work1ds, sendcounts, sdispls, & work1dr,recvcounts,sdispls,pelist) nullify(work1dr) nullify(work1ds) @@ -289,4 +295,6 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, enddo ! return - end + end subroutine sumfln_stochy + + end module sumfln_stochy_mod diff --git a/update_ca.f90 b/update_ca.f90 index 9bc38eae..10014d38 100644 --- a/update_ca.f90 +++ b/update_ca.f90 @@ -40,9 +40,9 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, logical :: ca_global, ca_sgs real :: Wpert(nxc,nyc),Cpert(nxc,nyc) -!SGS parameters: +!SGS parameters: integer(8) :: count, count_rate, count_max, count_trunc -integer(8) :: iscale = 10000000000 +integer(8) :: iscale = 10000000000 integer :: count5, count6 type(random_stat) :: rstate real :: dt, timescale, sigma, skew, kurt, acorr, gamma @@ -74,8 +74,8 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, if (.not. allocated(field_in))then allocate(field_in(nxc*nyc,1)) endif - if(.not. allocated(board_halo))then - allocate(board_halo(nxch,nych,1)) + if(.not. allocated(board_halo))then + allocate(board_halo(nxch,nych,1)) endif if(.not. allocated(Wpert_halo))then allocate(Wpert_halo(nxch,nych,1)) @@ -83,43 +83,43 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, if(.not. allocated(Cpert_halo))then allocate(Cpert_halo(nxch,nych,1)) endif - - if(ca_sgs == .true.)then + + if(ca_sgs)then if(kstep <= 1)then - + do j=1,nyc do i=1,nxc board(i,j,nf) = 0 lives(i,j,nf) = 0 enddo enddo - + endif if(kstep == 2)then !Initiate CA at kstep 2 as physics field is empty at 0 and 1. - + do j=1,nyc do i=1,nxc board(i,j,nf) = iini(i,j,nf) lives(i,j,nf) = ilives(i,j)*iini(i,j,nf) enddo enddo - + endif !Seed with new CA cells at each nseed step if(mod(kstep,nseed) == 0 .and. kstep >= 2)then - + do j=1,nyc do i=1,nxc - board(i,j,nf) = iini(i,j,nf) - lives(i,j,nf) = ilives(i,j)*iini(i,j,nf) + board(i,j,nf) = iini(i,j,nf) + lives(i,j,nf) = ilives(i,j)*iini(i,j,nf) enddo enddo - + endif - + if(kstep == 2)then spinup=nspinup else @@ -139,7 +139,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, endif - !Seed with new CA cells at each nseed step + !Seed with new CA cells at each nseed step if(mod(kstep,nseed) == 0)then do j=1,nyc @@ -156,18 +156,18 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, else spinup = 1 endif - + endif !sgs/global !Step 1 - Solve the stochastic gaussian skewed SGS equation in order to generate !perturbations to the grid mean model fields. -if (ca_sgs == .true.) then -!Compute the SGS and perturb the vertical velocity field: -!Read these values in from namelist, guided from LES data: +if (ca_sgs) then +!Compute the SGS and perturb the vertical velocity field: +!Read these values in from namelist, guided from LES data: -dt=900. +dt=900. timescale=21600.0 sigma=0.8 skew=0.8 @@ -176,7 +176,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, gamma=-log(acorr)/dt !calculate coeffients for SGS with auto-correlation and use -!these for predictor-correcting time stepping +!these for predictor-correcting time stepping E_sq2=2.0*gamma/3.0*(kurt-3.0/2.0*skew**2.)/(kurt-skew**2.+2.) E2=sqrt(E_sq2) g2D=skew*sigma*(gamma-E_sq2)/(2.0*E2) @@ -194,17 +194,17 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, do it=1,spinup -if (ca_sgs == .true.) then +if (ca_sgs) then - !Random seed for SGS + !Random seed for SGS noise1D1 = 0.0 noise1D2 = 0.0 - + call system_clock(count, count_rate, count_max) count_trunc = iscale*(count/iscale) count5 = count - count_trunc count6=count5+9827 - !broadcast to all tasks + !broadcast to all tasks call mp_bcst(count5) call mp_bcst(count6) @@ -213,7 +213,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, call random_setseed(count6,rstate) call random_gauss(noise1D2,rstate) - !Put on 2D: + !Put on 2D: do j=1,nyc do i=1,nxc NOISE_A(i,j)=noise1D1(i+(j-1)*nxc) @@ -260,7 +260,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, sgs2(i,j,nf)=sgs2(i,j,nf)-(flamx2*tmp1a+0.5*E2*g2D(i,j))*dt + (B2*NOISE_A(i,j)*sqrtdt) & + (g2D(i,j) + E2 * tmp1a )* NOISE_B(i,j)*sqrtdt - Wpert(i,j)=vertvelhigh(i,j)*(1.0 + sgs2(i,j,nf)) + Wpert(i,j)=vertvelhigh(i,j)*(1.0 + sgs2(i,j,nf)) enddo enddo @@ -268,9 +268,9 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, endif !ca sgs true - + !Step 2 - Initialize variables to 0 and extract the halo - + neighbours=0 birth=0 newlives=0 @@ -284,22 +284,22 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, Cpert_halo=0 field_in=0 -!The input to scalar_field_halo needs to be 1D. +!The input to scalar_field_halo needs to be 1D. !take the updated board fields and extract the halo -! in order to have updated values in the halo region. - - if(ca_global ==.true.)then - do j=1,nyc - do i=1,nxc - field_in(i+(j-1)*nxc,1)=board(i,j,nf) - enddo - enddo -!Step 3 - Extract the halo +! in order to have updated values in the halo region. + + if(ca_global)then + do j=1,nyc + do i=1,nxc + field_in(i+(j-1)*nxc,1)=board(i,j,nf) + enddo + enddo +!Step 3 - Extract the halo call atmosphere_scalar_field_halo(board_halo,halo,nxch,nych,k_in,field_in) endif - - if(ca_sgs==.true.)then + + if(ca_sgs)then field_in=0 do j=1,nyc do i=1,nxc @@ -307,7 +307,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, enddo enddo - call atmosphere_scalar_field_halo(Wpert_halo,halo,nxch,nych,k_in,field_in) + call atmosphere_scalar_field_halo(Wpert_halo,halo,nxch,nych,k_in,field_in) field_in=0 do j=1,nyc @@ -323,9 +323,9 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, !Step 4 - Compute the neighbourhood -if(ca_sgs == .true.)then !SGSmethod +if(ca_sgs)then !SGSmethod - !Count the number of neighbours where perturbed massflux is larger than + !Count the number of neighbours where perturbed massflux is larger than !a threshold @@ -386,7 +386,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, enddo enddo - + !CA stand alone method else !global @@ -403,7 +403,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, endif !sgs/global ! Step 5 - Check rules; the birth condition differs between SGS and GOL method -if(ca_sgs == .true.)then !SGS +if(ca_sgs)then !SGS if(nf==1)then do j=1,nyc @@ -483,7 +483,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, enddo enddo - if(ca_sgs == .true.)then + if(ca_sgs)then do j=1,nyc do i=1,nxc lives(i,j,nf)=lives(i,j,nf)+newcell(i,j)*ilives(i,j) @@ -510,7 +510,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, enddo !spinup !COARSE-GRAIN BACK TO NWP GRID - + inci=ncells incj=ncells sub=ncells-1 @@ -526,24 +526,24 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, lives_max=maxval(ilives) call mp_reduce_max(lives_max) - if(ca_sgs == .true.)then + if(ca_sgs)then CA(:,:) = (frac(:,:)/lives_max) else !global CA(:,:) = (frac(:,:)/real(nlives)) endif -if(nca_plumes == .true.) then +if(nca_plumes) then !COMPUTE NUMBER OF CLUSTERS (CONVECTIVE PLUMES) IN EACH CA-CELL !Note, at the moment we only use the count of the plumes found in a grid-cell -!In the future the routine "plumes" can also be used to give the size of +!In the future the routine "plumes" can also be used to give the size of !each individual plume for better coupling to the convection scheme. temp=0 do j=1,nyc do i=1,nxc if(lives(i,j,1) > 0)then - temp(i,j)=1 + temp(i,j)=1 endif enddo enddo @@ -555,7 +555,7 @@ subroutine update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini, if (.not. allocated(L))then allocate(L(kend)) endif - + ca_plumes(:,:)=0 inci=ncells incj=ncells