From ec86a8ec218aba1a7b77ad5aed46e45925da93ec Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 3 Dec 2019 21:58:18 +0000 Subject: [PATCH 01/10] CA global updates --- cellular_automata_global.F90 | 484 ++++++++++++++ ..._automata.f90 => cellular_automata_sgs.F90 | 429 ++++++------ makefile | 7 +- update_ca.F90 | 608 ++++++++++++++++++ update_ca.f90 | 586 ----------------- 5 files changed, 1326 insertions(+), 788 deletions(-) create mode 100644 cellular_automata_global.F90 rename cellular_automata.f90 => cellular_automata_sgs.F90 (52%) create mode 100644 update_ca.F90 delete mode 100644 update_ca.f90 diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 new file mode 100644 index 0000000..1f67a76 --- /dev/null +++ b/cellular_automata_global.F90 @@ -0,0 +1,484 @@ +subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & + nca,ncells,nlives,nfracseed,nseed,nthresh,ca_global,ca_sgs,iseed_ca, & + ca_smooth,nsmooth,nspinup,blocksize) + +use machine +use update_ca, only: update_cells_sgs, update_cells_global +#ifdef STOCHY_UNIT_TEST +use standalone_stochy_module, only: GFS_Coupling_type, GFS_diag_type, GFS_statein_type +use atmosphere_stub_mod, only: atmosphere_resolution, atmosphere_domain, & + atmosphere_scalar_field_halo, atmosphere_control_data +#else +use GFS_typedefs, only: GFS_Coupling_type, GFS_diag_type, GFS_statein_type +use atmosphere_mod, only: atmosphere_resolution, atmosphere_domain, & + atmosphere_scalar_field_halo, atmosphere_control_data +#endif +use mersenne_twister, only: random_setseed,random_gauss,random_stat,random_number +use mpp_domains_mod, only: domain2D +use block_control_mod, only: block_control_type, define_blocks_packed +use fv_mp_mod, only : mp_reduce_sum,mp_bcst,mp_reduce_max,mp_reduce_min,is_master + + +implicit none + +!L.Bengtsson, 2017-06 + +!This program evolves a cellular automaton uniform over the globe given +!the flag ca_global + +integer,intent(in) :: kstep,ncells,nca,nlives,nseed,iseed_ca,nspinup,nsmooth +real,intent(in) :: nfracseed,nthresh +logical,intent(in) :: ca_global, ca_sgs, ca_smooth +integer, intent(in) :: nblks,nlev,blocksize +type(GFS_coupling_type),intent(inout) :: Coupling(nblks) +type(GFS_diag_type),intent(inout) :: Diag(nblks) +type(GFS_statein_type),intent(in) :: Statein(nblks) +type(block_control_type) :: Atm_block +type(random_stat) :: rstate +integer :: nlon, nlat, isize,jsize,nf,nn +integer :: inci, incj, nxc, nyc, nxch, nych +integer :: halo, k_in, i, j, k, iec, jec, isc, jsc +integer :: seed, ierr7,blk, ix, iix, count4,ih,jh +integer :: blocksz,levs +integer(8) :: count, count_rate, count_max, count_trunc +integer(8) :: iscale = 10000000000 +integer, allocatable :: iini(:,:,:),ilives(:,:),iini_g(:,:,:),ilives_g(:,:),ca_plumes(:,:) +real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:),Detfield(:,:,:) +real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:) +real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:),tmp(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),condition(:,:),rho(:,:),cape(:,:) +real(kind=kind_phys), allocatable :: noise1D(:),vertvelhigh(:,:),noise(:,:,:) +real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,lambda +real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt +logical,save :: block_message=.true. +logical :: nca_plumes + +!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. +!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. +!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. + +halo=1 +k_in=1 + +nca_plumes = .false. +!---------------------------------------------------------------------------- +! Get information about the compute domain, allocate fields on this +! domain + +! WRITE(*,*)'Entering cellular automata calculations' + +! Some security checks for namelist combinations: + if(nca > 5)then + write(0,*)'Namelist option nca cannot be larger than 5 - exiting' + stop + endif + +! if(ca_global == .true. .and. ca_sgs == .true.)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 +! write(0,*)'Currently ca_smooth does not work with ca_sgs - exiting' +! stop +! endif + + 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. + + inci=ncells + incj=ncells + + nxc=nlon*ncells + nyc=nlat*ncells + + nxch=nxc+2*halo + nych=nyc+2*halo + + + !Allocate fields: + + allocate(cloud(nlon,nlat)) + allocate(omega(nlon,nlat,29)) + allocate(pressure(nlon,nlat,29)) + allocate(humidity(nlon,nlat)) + allocate(dp(nlon,nlat,29)) + allocate(rho(nlon,nlat)) + allocate(surfp(nlon,nlat)) + allocate(vertvelmean(nlon,nlat)) + allocate(vertvelsum(nlon,nlat)) + allocate(field_in(nlon*nlat,1)) + allocate(field_out(isize,jsize,1)) + allocate(field_smooth(nlon,nlat)) + allocate(iini(nxc,nyc,nca)) + allocate(ilives(nxc,nyc)) + allocate(iini_g(nxc,nyc,nca)) + allocate(ilives_g(nxc,nyc)) + allocate(vertvelhigh(nxc,nyc)) + allocate(condition(nxc,nyc)) + allocate(cape(nlon,nlat)) + allocate(Detfield(nlon,nlat,nca)) + allocate(CA(nlon,nlat)) + allocate(ca_plumes(nlon,nlat)) + allocate(CA1(nlon,nlat)) + allocate(CA2(nlon,nlat)) + allocate(CA3(nlon,nlat)) + allocate(tmp(nlon,nlat)) + allocate(noise(nxc,nyc,nca)) + allocate(noise1D(nxc*nyc)) + + !Initialize: + Detfield(:,:,:)=0. + vertvelmean(:,:) =0. + vertvelsum(:,:)=0. + cloud(:,:)=0. + humidity(:,:)=0. + condition(:,:)=0. + cape(:,:)=0. + vertvelhigh(:,:)=0. + noise(:,:,:) = 0.0 + noise1D(:) = 0.0 + iini(:,:,:) = 0 + ilives(:,:) = 0 + iini_g(:,:,:) = 0 + ilives_g(:,:) = 0 + CA1(:,:) = 0.0 + CA2(:,:) = 0.0 + CA3(:,:) = 0.0 + ca_plumes(:,:) = 0 + +!Put the blocks of model fields into a 2d array + levs=nlev + blocksz=blocksize + + call atmosphere_control_data(isc, iec, jsc, jec, levs) + call define_blocks_packed('cellular_automata', Atm_block, isc, iec, jsc, jec, levs, & + blocksz, block_message) + + isc = Atm_block%isc + iec = Atm_block%iec + jsc = Atm_block%jsc + 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) + surfp(i,j) = Statein(blk)%pgr(ix) + humidity(i,j)=Statein(blk)%qgrs(ix,13,1) !about 850 hpa + do k = 1,29 !Lower troposphere: level 29 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 + enddo + enddo + +!Compute layer averaged vertical velocity (Pa/s) + vertvelsum=0. + vertvelmean=0. + do j=1,nlat + do i =1,nlon + dp(i,j,1)=(surfp(i,j)-pressure(i,j,1)) + do k=2,29 + dp(i,j,k)=(pressure(i,j,k-1)-pressure(i,j,k)) + enddo + count1=0. + do k=1,29 + count1=count1+1. + vertvelsum(i,j)=vertvelsum(i,j)+(omega(i,j,k)*dp(i,j,k)) + enddo + enddo + enddo + + do j=1,nlat + do i=1,nlon + vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,29)) + enddo + enddo + +!Generate random number, following stochastic physics code: +do nf=1,nca + + if (iseed_ca == 0) then + ! generate a random seed from system clock and ens member number + call system_clock(count, count_rate, count_max) + ! iseed is elapsed time since unix epoch began (secs) + ! truncate to 4 byte integer + count_trunc = iscale*(count/iscale) + count4 = count - count_trunc + nf*201 + else + ! don't rely on compiler to truncate integer(8) to integer(4) on + ! overflow, do wrap around explicitly. + count4 = mod(iseed_ca + 2147483648, 4294967296) - 2147483648 + nf*201 + endif + + !Set seed (to be the same) on all tasks. Save random state. + call random_setseed(count4,rstate) + call random_number(noise1D,rstate) + !Put on 2D: + do j=1,nyc + do i=1,nxc + noise(i,j,nf)=noise1D(i+(j-1)*nxc) + enddo + enddo + + +!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 + iini_g(i,j,nf)=1 + else + iini_g(i,j,nf)=0 + endif + enddo + 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.) + + +do nf=1,nca !update each ca + + do j = 1,nyc + do i = 1,nxc + ilives_g(i,j)=int(real(nlives)*1.5*noise(i,j,nf)) + enddo + enddo + + + +!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 + + + CA(:,:)=0. + + call update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini_g,ilives_g, & + nlives, ncells, nfracseed, nseed,nthresh, ca_global, & + ca_sgs,nspinup, condition, vertvelhigh,nf,nca_plumes) + + + +if (ca_smooth) then +do nn=1,nsmooth !number of itterations for the smoothing. + +field_in=0. + +!get halo +do j=1,nlat + do i=1,nlon + field_in(i+(j-1)*nlon,1)=CA(i,j) + enddo +enddo + +field_out=0. + +call atmosphere_scalar_field_halo(field_out,halo,isize,jsize,k_in,field_in) + +do j=1,nlat + 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)+ & + 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)+& + 2.0*field_out(ih+1,jh-1,1))/20. + enddo +enddo + +do j=1,nlat + do i=1,nlon + CA(i,j)=field_smooth(i,j) + enddo +enddo + +enddo !nn +endif + + if(nf==1)then + CA1(:,:)=CA(:,:) + elseif(nf==2)then + CA2(:,:)=CA(:,:) + else + CA3(:,:)=CA(:,:) + endif + + +enddo !nca loop + + +!!!!Post processing, should be made into a separate subroutine + +!!!!Construct linear combinations of CA1, CA2 and CA3 + + +!Use min-max method to normalize range +!!!CA1 +Detmax(1)=maxval(CA1) +call mp_reduce_max(Detmax(1)) +Detmin(1)=minval(CA1) +call mp_reduce_min(Detmin(1)) + +do j=1,nlat + do i=1,nlon + CA1(i,j) = ((CA1(i,j) - Detmin(1))/(Detmax(1)-Detmin(1))) + enddo +enddo + +CAmean=0. +psum=0. +csum=0. +do j=1,nlat + do i=1,nlon + psum=psum+(CA1(i,j)) + csum=csum+1 + enddo +enddo + +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) + +CAmean=psum/csum + +do j=1,nlat + do i=1,nlon + CA1(i,j)=(CA1(i,j)-CAmean)+1.0 !Can we compute the global median? + enddo +enddo + + + +!!!CA2 +Detmax(2)=maxval(CA2) +call mp_reduce_max(Detmax(2)) +Detmin(2)=minval(CA2) +call mp_reduce_min(Detmin(2)) + +do j=1,nlat + do i=1,nlon + CA2(i,j) = ((CA2(i,j) - Detmin(2))/(Detmax(2)-Detmin(2))) + enddo +enddo + +CAmean=0. +psum=0. +csum=0. +do j=1,nlat + do i=1,nlon + psum=psum+(CA2(i,j)) + csum=csum+1 + enddo +enddo + +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) + +CAmean=psum/csum + +do j=1,nlat + do i=1,nlon + CA2(i,j)=(CA2(i,j)-CAmean)+1.0 !Can we compute the global median? + enddo +enddo + + +!!!CA3 +Detmax(3)=maxval(CA3) +call mp_reduce_max(Detmax(3)) +Detmin(3)=minval(CA3) +call mp_reduce_min(Detmin(3)) + +do j=1,nlat + do i=1,nlon + CA3(i,j) = ((CA3(i,j) - Detmin(3))/(Detmax(3)-Detmin(3))) + enddo +enddo + +CAmean=0. +psum=0. +csum=0. +do j=1,nlat + do i=1,nlon + psum=psum+(CA3(i,j)) + csum=csum+1 + enddo +enddo + +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) + +CAmean=psum/csum + +do j=1,nlat + do i=1,nlon + CA3(i,j)=(CA3(i,j)-CAmean)+1.0 !Can we compute the global median? + enddo +enddo + + + +!Put back into blocks 1D array to be passed to physics +!or diagnostics output +if(kstep < 1)then +CA1(:,:)=1. +CA2(:,:)=1. +CA3(:,:)=1. +endif + + + 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 + Diag(blk)%ca1(ix)=CA1(i,j) + Diag(blk)%ca2(ix)=CA2(i,j) + Diag(blk)%ca3(ix)=CA3(i,j) + Coupling(blk)%ca1(ix)=CA1(i,j) + Coupling(blk)%ca2(ix)=CA2(i,j) + Coupling(blk)%ca3(ix)=CA3(i,j) + enddo + enddo + + deallocate(omega) + deallocate(pressure) + deallocate(humidity) + deallocate(dp) + deallocate(cape) + deallocate(rho) + deallocate(surfp) + deallocate(vertvelmean) + deallocate(vertvelsum) + deallocate(field_in) + deallocate(field_out) + deallocate(field_smooth) + deallocate(iini) + deallocate(ilives) + deallocate(condition) + deallocate(Detfield) + deallocate(CA) + deallocate(ca_plumes) + deallocate(CA1) + deallocate(CA2) + deallocate(CA3) + deallocate(noise) + deallocate(noise1D) + +end subroutine cellular_automata_global diff --git a/cellular_automata.f90 b/cellular_automata_sgs.F90 similarity index 52% rename from cellular_automata.f90 rename to cellular_automata_sgs.F90 index b8da0c2..b0a7822 100644 --- a/cellular_automata.f90 +++ b/cellular_automata_sgs.F90 @@ -1,16 +1,22 @@ -subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & +subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & nca,ncells,nlives,nfracseed,nseed,nthresh,ca_global,ca_sgs,iseed_ca, & ca_smooth,nspinup,blocksize) use machine -use update_ca, only: update_cells +use update_ca, only: update_cells_sgs, update_cells_global +#ifdef STOCHY_UNIT_TEST +use standalone_stochy_module, only: GFS_Coupling_type, GFS_diag_type, GFS_statein_type +use atmosphere_stub_mod, only: atmosphere_resolution, atmosphere_domain, & + atmosphere_scalar_field_halo, atmosphere_control_data +#else +use GFS_typedefs, only: GFS_Coupling_type, GFS_diag_type, GFS_statein_type use atmosphere_mod, only: atmosphere_resolution, atmosphere_domain, & atmosphere_scalar_field_halo, atmosphere_control_data +#endif use mersenne_twister, only: random_setseed,random_gauss,random_stat,random_number -use GFS_typedefs, only: GFS_Coupling_type, GFS_diag_type, GFS_statein_type use mpp_domains_mod, only: domain2D use block_control_mod, only: block_control_type, define_blocks_packed -use fv_mp_mod, only : mp_reduce_sum,mp_bcst,mp_reduce_max,is_master +use fv_mp_mod, only : mp_reduce_sum,mp_bcst,mp_reduce_max,mp_reduce_min,is_master implicit none @@ -18,20 +24,18 @@ 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 -!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 @@ -42,44 +46,37 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & type(GFS_statein_type),intent(in) :: Statein(nblks) type(block_control_type) :: Atm_block type(random_stat) :: rstate -integer :: nlon, nlat, isize,jsize,nf,k350,k850 +integer :: nlon, nlat, isize,jsize,nf,nn integer :: inci, incj, nxc, nyc, nxch, nych integer :: halo, k_in, i, j, k, iec, jec, isc, jsc integer :: seed, ierr7,blk, ix, iix, count4,ih,jh -integer :: blocksz,levs,ra,rb,rc +integer :: blocksz,levs integer(8) :: count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 -integer, allocatable :: iini(:,:,:),ilives(:,:),ca_plumes(:,:) +integer, allocatable :: iini(:,:,:),ilives(:,:,:),iini_g(:,:,:),ilives_g(:,:),ca_plumes(:,:) real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:),Detfield(:,:,:) -real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:) -real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:) -real(kind=kind_phys), allocatable :: CA(:,:),CAstore(:,:),CAavg(:,:),condition(:,:),rho(:,:),cape(:,:) -real(kind=kind_phys), allocatable :: CA_DEEP(:,:),CA_TURB(:,:),CA_RAD(:,:),CA_MICRO(:,:),CA_SHAL(:,:) +real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:),uwind(:,:) +real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:),shalp(:,:),gamt(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),rho(:,:),cape(:,:) +real(kind=kind_phys), allocatable :: CA_DEEP(:,:),CA_TURB(:,:),CA_SHAL(:,:) real(kind=kind_phys), allocatable :: noise1D(:),vertvelhigh(:,:),noise(:,:,:) -real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,alpha -real(kind=kind_phys) :: Detmax(nca),Detmean(nca),phi,stdev,delt +real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,lambda +real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt logical,save :: block_message=.true. logical :: nca_plumes !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. -k350 = 29 -k850 = 13 -alpha = 1.5 -ra = 201 -rb = 2147483648 -rc = 4294967296 - halo=1 k_in=1 @@ -88,47 +85,33 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & ! Get information about the compute domain, allocate fields on this ! domain -! WRITE(*,*)'Entering cellular automata calculations' - ! Some security checks for namelist combinations: if(nca > 5)then write(0,*)'Namelist option nca cannot be larger than 5 - exiting' stop endif - 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 .and. ca_smooth)then - write(0,*)'Currently ca_smooth does not work with ca_sgs - exiting' - stop - endif - 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. inci=ncells incj=ncells - + nxc=nlon*ncells nyc=nlat*ncells - + nxch=nxc+2*halo nych=nyc+2*halo - !Allocate fields: allocate(cloud(nlon,nlat)) - allocate(omega(nlon,nlat,k350)) - allocate(pressure(nlon,nlat,k350)) + allocate(omega(nlon,nlat,29)) + allocate(pressure(nlon,nlat,29)) allocate(humidity(nlon,nlat)) - allocate(dp(nlon,nlat,k350)) + allocate(uwind(nlon,nlat)) + allocate(dp(nlon,nlat,29)) allocate(rho(nlon,nlat)) allocate(surfp(nlon,nlat)) allocate(vertvelmean(nlon,nlat)) @@ -137,44 +120,42 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & allocate(field_out(isize,jsize,1)) allocate(field_smooth(nlon,nlat)) allocate(iini(nxc,nyc,nca)) - allocate(ilives(nxc,nyc)) + allocate(ilives(nxc,nyc,nca)) + allocate(iini_g(nxc,nyc,nca)) + allocate(ilives_g(nxc,nyc)) allocate(vertvelhigh(nxc,nyc)) allocate(condition(nxc,nyc)) allocate(cape(nlon,nlat)) + allocate(shalp(nlon,nlat)) + allocate(gamt(nlon,nlat)) allocate(Detfield(nlon,nlat,nca)) allocate(CA(nlon,nlat)) allocate(ca_plumes(nlon,nlat)) - allocate(CAstore(nlon,nlat)) - allocate(CAavg(nlon,nlat)) allocate(CA_TURB(nlon,nlat)) - allocate(CA_RAD(nlon,nlat)) - allocate(CA_DEEP(nlon,nlat)) - allocate(CA_MICRO(nlon,nlat)) + allocate(CA_DEEP(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. + uwind(:,:) = 0. condition(:,:)=0. cape(:,:)=0. vertvelhigh(:,:)=0. + ca_plumes(:,:) = 0 noise(:,:,:) = 0.0 noise1D(:) = 0.0 iini(:,:,:) = 0 - ilives(:,:) = 0 - CA(:,:) = 0.0 - CAavg(:,:) = 0.0 - ca_plumes(:,:) = 0 - CA_TURB(:,:) = 0.0 - CA_RAD(:,:) = 0.0 - CA_DEEP(:,:) = 0.0 - CA_MICRO(:,:) = 0.0 - CA_SHAL(:,:) = 0.0 + ilives(:,:,:) = 0 + iini_g(:,:,:) = 0 + ilives_g(:,:) = 0 + Detmax(:)=0. + Detmin(:)=0. !Put the blocks of model fields into a 2d array levs=nlev @@ -187,33 +168,40 @@ 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) + uwind(i,j) = Statein(blk)%ugrs(ix,30) + cape(i,j) = Coupling(blk)%cape(ix) + shalp(i,j) = Coupling(blk)%ca_shal(ix) + gamt(i,j) = Coupling(blk)%ca_turb(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 + humidity(i,j)=Statein(blk)%qgrs(ix,13,1) !about 850 hpa + do k = 1,29 !Lower troposphere: level 29 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 enddo enddo + CA_TURB(:,:) = 0.0 + CA_DEEP(:,:) = 0.0 + CA_SHAL(:,:) = 0.0 + !Compute layer averaged vertical velocity (Pa/s) vertvelsum=0. vertvelmean=0. do j=1,nlat do i =1,nlon dp(i,j,1)=(surfp(i,j)-pressure(i,j,1)) - do k=2,k350 + do k=2,29 dp(i,j,k)=(pressure(i,j,k-1)-pressure(i,j,k)) enddo count1=0. - do k=1,k350 + do k=1,29 count1=count1+1. vertvelsum(i,j)=vertvelsum(i,j)+(omega(i,j,k)*dp(i,j,k)) enddo @@ -222,7 +210,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon - vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,k350)) + vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,29)) enddo enddo @@ -235,26 +223,25 @@ 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*201 else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. - count4 = mod(iseed_ca + rb, rc) - rb + nf*ra + count4 = mod(iseed_ca + 2147483648, 4294967296) - 2147483648 + nf*201 endif - !Set seed (to be the same) on all tasks. Save random state. - call random_setseed(count4,rstate) - call random_number(noise1D,rstate) + call random_setseed(count4) + call random_number(noise1D) !Put on 2D: do j=1,nyc do i=1,nxc - noise(i,j,nf)=noise1D(i+(j-1)*nxc) + noise(i,j,nf)=noise1D(i+(j-1)*nxc) enddo enddo !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 @@ -264,25 +251,23 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo 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.) -CAstore = 0. do nf=1,nca !update each ca - - if(ca_sgs)then - - if(nf==1)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) !cape(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -294,8 +279,8 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo do j = 1,nyc - do i = 1,nxc - ilives(i,j)=int(real(nlives)*alpha*noise(i,j,nf)) + do i = 1,nxc + ilives(i,j,nf)=nlives enddo enddo @@ -304,7 +289,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=humidity(inci/ncells,incj/ncells) + condition(i,j)=shalp(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -317,7 +302,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,nf)=nlives !int(real(nlives)*1.5*noise(i,j,nf)) enddo enddo @@ -326,7 +311,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) + condition(i,j)=gamt(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -339,7 +324,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,nf)=nlives enddo enddo @@ -348,7 +333,7 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) + condition(i,j)=vertvelmean(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -361,13 +346,13 @@ 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,nf)=int(real(nlives)*1.5*noise(i,j,nf)) enddo enddo else inci=ncells - incj=ncells + incj=ncells do j=1,nyc do i=1,nxc condition(i,j)=cape(inci/ncells,incj/ncells) @@ -383,14 +368,14 @@ 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,nf)=int(real(nlives)*1.5*noise(i,j,nf)) enddo enddo endif !nf - + !Vertical velocity has its own variable in order to condition on combination !of "condition" and vertical velocity. @@ -409,144 +394,191 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo - - else !ca_global - - do j = 1,nyc - do i = 1,nxc - ilives(i,j)=int(real(nlives)*alpha*noise(i,j,nf)) - enddo - enddo - - 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 - - call update_cells(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini,ilives, & +!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_sgs(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)then - CAstore(:,:) = CAstore(:,:) + CA(:,:) - elseif(ca_sgs)then if(nf==1)then - CA_DEEP(:,:)=CA(:,:) + CA_DEEP(:,:)=CA(:,:) elseif(nf==2)then - CA_TURB(:,:)=CA(:,:) - elseif(nf==3)then - CA_SHAL(:,:)=CA(:,:) - elseif(nf==4)then - CA_RAD(:,:)=CA(:,:) + CA_SHAL(:,:)=CA(:,:) else - CA_MICRO(:,:)=CA(:,:) + CA_TURB(:,:)=CA(:,:) endif - else - write(*,*)'Either sgs or global needs to be selected' - endif + enddo !nf (nca) - if(ca_global)then - CAavg = CAstore / real(nca) - endif +!!Post-processesing - could be made into a separate sub-routine -!smooth CA field +!Deep convection ==================================== -if (ca_smooth .and. ca_global) then -field_in=0. +if(kstep > 1)then -!get halo +!Box-Cox method to make output range less skewed +lambda=0.38 do j=1,nlat do i=1,nlon - field_in(i+(j-1)*nlon,1)=CAavg(i,j) + if(CA_DEEP(i,j).NE.0.)then + CA_DEEP(i,j)=((CA_DEEP(i,j)**lambda)-1.0)/lambda + endif enddo enddo -field_out=0. +!Use min-max method to normalize range +Detmax(1)=maxval(CA_DEEP,CA_DEEP.NE.0.) +call mp_reduce_max(Detmax(1)) +Detmin(1)=minval(CA_DEEP,CA_DEEP.NE.0.) +call mp_reduce_min(Detmin(1)) -call atmosphere_scalar_field_halo(field_out,halo,isize,jsize,k_in,field_in) do j=1,nlat 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)+ & - 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)+& - 2.0*field_out(ih+1,jh-1,1))/20. + if(CA_DEEP(i,j).NE.0.)then + CA_DEEP(i,j) =(CA_DEEP(i,j) - Detmin(1))/(Detmax(1)-Detmin(1)) + endif enddo enddo +!Compute the mean of the new range and subtract +CAmean=0. +psum=0. +csum=0. do j=1,nlat do i=1,nlon - CAavg(i,j)=field_smooth(i,j) + if(CA_DEEP(i,j).NE.0.)then + psum=psum+(CA_DEEP(i,j)) + csum=csum+1 + endif enddo enddo -endif !smooth +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) -!In case of ca_global give data zero mean and unit standard deviation +CAmean=psum/csum -!if(ca_global == .true.)then +do j=1,nlat + do i=1,nlon + if(CA_DEEP(i,j).NE.0.)then + CA_DEEP(i,j)=(CA_DEEP(i,j)-CAmean) + endif + enddo +enddo -!CAmean=0. -!psum=0. -!csum=0. +Detmin(1) = minval(CA_DEEP,CA_DEEP.NE.0) +call mp_reduce_min(Detmin(1)) + +!Set range outside positive vertical mean vertical velocity to 0. !do j=1,nlat ! do i=1,nlon -! psum=psum+CAavg(i,j) -! csum=csum+1 +! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward +! CA_DEEP(i,j)=0. +! endif ! enddo !enddo -!call mp_reduce_sum(psum) -!call mp_reduce_sum(csum) -!CAmean=psum/csum -!CAstdv=0. -!sq_diff = 0. -!do j=1,nlat -! do i=1,nlon -! sq_diff = sq_diff + (CAavg(i,j)-CAmean)**2.0 -! enddo -!enddo +!Shallow convection ============================================================ +!Box-Cox method to make range less skewed + +lambda=0.38 +do j=1,nlat + do i=1,nlon + if(CA_SHAL(i,j).NE.0.)then + CA_SHAL(i,j)=((CA_SHAL(i,j)**lambda)-1.0)/lambda + endif + enddo +enddo + +!Use min-max method to normalize range +Detmax(2)=maxval(CA_SHAL,CA_SHAL.NE.0) +call mp_reduce_max(Detmax(2)) +Detmin(2)=minval(CA_SHAL,CA_SHAL.NE.0) +call mp_reduce_min(Detmin(2)) + +do j=1,nlat + do i=1,nlon + if(CA_SHAL(i,j).NE.0.)then + CA_SHAL(i,j)=(CA_SHAL(i,j) - Detmin(2))/(Detmax(2)-Detmin(2)) + endif + enddo +enddo + +!Compute the mean of the new range and subtract +CAmean=0. +psum=0. +csum=0. +do j=1,nlat + do i=1,nlon + if(CA_SHAL(i,j).NE.0.)then + psum=psum+(CA_SHAL(i,j)) + csum=csum+1 + endif + enddo +enddo -!call mp_reduce_sum(sq_diff) +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) -!CAstdv = sqrt( sq_diff / (csum-1.0) ) +CAmean=psum/csum +do j=1,nlat + do i=1,nlon + if(CA_SHAL(i,j).NE.0.)then + CA_SHAL(i,j)=(CA_SHAL(i,j)-CAmean) + endif + enddo +enddo + +!Set range outside positive vertical mean vertical velocity to 0. !do j=1,nlat ! do i=1,nlon -! CAavg(i,j)=(CAavg(i,j)-CAmean)/CAstdv +! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward ! +! CA_SHAL(i,j)=0. +! endif ! enddo !enddo -!endif -!Set the range for the nca individual ca_sgs patterns: -if(ca_sgs)then - -Detmax(1)=maxval(CA_DEEP(:,:)) -call mp_reduce_max(Detmax(1)) +!Turbulence ============================================================================= +!Box-Cox method to make the range less skewed +lambda=0.38 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 + if(CA_TURB(i,j).NE.0.)then + CA_TURB(i,j)=((CA_TURB(i,j)**lambda)-1.0)/lambda endif enddo enddo +!Use min-max method to normalize range +Detmax(3)=maxval(CA_TURB,CA_TURB.NE.0) +call mp_reduce_max(Detmax(3)) +Detmin(3)=minval(CA_TURB,CA_TURB.NE.0) +call mp_reduce_min(Detmin(3)) + +do j=1,nlat + do i=1,nlon + if(CA_TURB(i,j).NE.0.)then + CA_TURB(i,j)=(CA_TURB(i,j) - Detmin(3))/(Detmax(3)-Detmin(3)) + endif + enddo +enddo + +!Compute the mean of the new range and subtract CAmean=0. psum=0. csum=0. do j=1,nlat do i=1,nlon - if(CA_DEEP(i,j)>0.)then - psum=psum+(CA_DEEP(i,j)) + if(CA_TURB(i,j).NE.0.)then + psum=psum+(CA_TURB(i,j)) csum=csum+1 endif enddo @@ -559,17 +591,25 @@ 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? + if(CA_TURB(i,j).NE.0.)then + CA_TURB(i,j)=(CA_TURB(i,j)-CAmean) endif enddo enddo -!!! +!Set range outside positive vertical mean vertical velocity to 0. +!do j=1,nlat +! do i=1,nlon +! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward ! +! CA_TURB(i,j)=0. +! endif +! enddo +!enddo +endif !kstep >1 -!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 @@ -578,35 +618,30 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo -endif - !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 j = Atm_block%index(blk)%jj(ix) - jsc + 1 - Diag(blk)%ca_out(ix)=CAavg(i,j) Diag(blk)%ca_deep(ix)=CA_DEEP(i,j) Diag(blk)%ca_turb(ix)=CA_TURB(i,j) Diag(blk)%ca_shal(ix)=CA_SHAL(i,j) - Diag(blk)%ca_rad(ix)=CA_RAD(i,j) - Diag(blk)%ca_micro(ix)=CA_MICRO(i,j) - Coupling(blk)%ca_out(ix)=CAavg(i,j) !Cellular Automata Coupling(blk)%ca_deep(ix)=CA_DEEP(i,j) Coupling(blk)%ca_turb(ix)=CA_TURB(i,j) Coupling(blk)%ca_shal(ix)=CA_SHAL(i,j) - Coupling(blk)%ca_rad(ix)=CA_RAD(i,j) - Coupling(blk)%ca_micro(ix)=CA_MICRO(i,j) enddo enddo + deallocate(omega) deallocate(pressure) deallocate(humidity) deallocate(dp) deallocate(cape) + deallocate(shalp) + deallocate(gamt) deallocate(rho) deallocate(surfp) deallocate(vertvelmean) @@ -620,14 +655,10 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & deallocate(Detfield) deallocate(CA) deallocate(ca_plumes) - deallocate(CAstore) - deallocate(CAavg) deallocate(CA_TURB) deallocate(CA_DEEP) - deallocate(CA_MICRO) deallocate(CA_SHAL) - deallocate(CA_RAD) deallocate(noise) deallocate(noise1D) -end subroutine cellular_automata +end subroutine cellular_automata_sgs diff --git a/makefile b/makefile index c863929..40060ec 100644 --- a/makefile +++ b/makefile @@ -22,8 +22,6 @@ FFLAGS += -I./include -I../FV3/gfsphysics/ -I../FV3/atmos_cubed_sphere -I$(FMS SRCS_F = SRCS_f90 = \ - ./cellular_automata.f90 \ - ./update_ca.f90 \ ./plumes.f90 SRCS_f = \ @@ -55,7 +53,10 @@ SRCS_F90 = \ ./stochy_patterngenerator.F90 \ ./stochy_data_mod.F90 \ ./get_stochy_pattern.F90 \ - ./initialize_spectral_mod.F90 + ./initialize_spectral_mod.F90 \ + ./cellular_automata_global.F90 \ + ./cellular_automata_sgs.F90 \ + ./update_ca.F90 SRCS_c = diff --git a/update_ca.F90 b/update_ca.F90 new file mode 100644 index 0000000..106434a --- /dev/null +++ b/update_ca.F90 @@ -0,0 +1,608 @@ +module update_ca + +#ifdef STOCHY_UNIT_TEST +use atmosphere_stub_mod, only: atmosphere_scalar_field_halo +#else +use atmosphere_mod, only: atmosphere_scalar_field_halo +#endif +use mersenne_twister, only: random_setseed,random_gauss,random_stat,random_number +use fv_mp_mod, only : mp_reduce_sum,mp_bcst,mp_reduce_min,mp_reduce_max + +implicit none + +!L. Bengtsson 2017-06 +!Evolve the cellular automata in time + + +contains + +subroutine update_cells_sgs(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) + +implicit none + +integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca +integer, intent(in) :: iini(nxc,nyc,nca) +integer, intent(inout) :: ilives(nxc,nyc,nca) +real, intent(in) :: condition(nxc,nyc),vertvelhigh(nxc,nyc) +real, intent(out) :: CA(nlon,nlat) +integer,intent(out) :: ca_plumes(nlon,nlat) +integer, intent(in) :: nlives, ncells, nseed, nspinup, nf +real, intent(in) :: nfracseed, nthresh +logical,intent(in) :: nca_plumes +real, dimension(nlon,nlat) :: frac +integer, dimension(nlon,nlat) :: maxlives +integer,allocatable,save :: board(:,:,:), lives(:,:,:) +integer,allocatable :: V(:),L(:) +integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize +integer :: haloh, ih, jh,lives_max,kend +real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh +real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) +integer, dimension(nxc,nyc) :: neighbours, birth, newlives,thresh,maxliveshigh +integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp +integer, dimension(ncells,ncells) :: onegrid +real,dimension(nxc,nyc) :: livesout +logical :: ca_global, ca_sgs +real :: Wpert(nxc,nyc),Cpert(nxc,nyc) + +!SGS parameters: +integer(8) :: count, count_rate, count_max, count_trunc +integer(8) :: iscale = 10000000000 +integer :: count5, count6 +type(random_stat) :: rstate +real :: dt, timescale, sigma, skew, kurt, acorr, gamma +real :: E_sq2, E2, g2, B_sq2, B2, sqrtdt,flamx2, tmp1, tmp1a +real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D +real, dimension(nxc*nyc) :: noise1D2, noise1D1 +real, allocatable, save :: sgs1(:,:,:),sgs2(:,:,:) +integer, dimension(nxch,nych,1) :: M_halo + + +!------------------------------------------------------------------------------------------------- +halo=1 +isize=nlon+2*halo +jsize=nlat+2*halo +k_in=1 + + if (.not. allocated(board))then + allocate(board(nxc,nyc,nca)) + endif + if (.not. allocated(lives))then + allocate(lives(nxc,nyc,nca)) + endif + if (.not. allocated(sgs1))then + allocate(sgs1(nxc,nyc,nca)) + endif + if (.not. allocated(sgs2))then + allocate(sgs2(nxc,nyc,nca)) + endif + 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)) + endif + if(.not. allocated(Wpert_halo))then + allocate(Wpert_halo(nxch,nych,1)) + endif + if(.not. allocated(Cpert_halo))then + allocate(Cpert_halo(nxch,nych,1)) + endif + +! Some initial white noise generation: (could be 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 + !call mp_bcst(count5) + !call mp_bcst(count6) + + call random_setseed(count5) + call random_number(noise1D1) + + call random_setseed(count6) + call random_number(noise1D2) + + !Put on 2D: +! do j=1,nyc +! do i=1,nxc +! NOISE_A(i,j)=noise1D1(i+(j-1)*nxc)-0.5 +! NOISE_B(i,j)=noise1D2(i+(j-1)*nxc) +! enddo +! enddo + +!Wpert=0. +!Cpert=0. + +! do j=1,nyc +! do i=1,nxc +! if(vertvelhigh(i,j) < 0. )then !upward motion +! Wpert(i,j)=vertvelhigh(i,j)*(1.0 + NOISE_A(i,j)) +! endif +! enddo +! enddo + +!endif + + +! do j=1,nyc +! do i=1,nxc +! if(vertvelhigh(i,j) < 0.)then !upward vertical motion +! thresh(i,j)=vertvelhigh(i,j) !condition(i,j) +! endif +! enddo +! enddo + + 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,nf)*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,nf)*iini(i,j,nf) + enddo + enddo + + endif + + if(kstep == 2)then + spinup=nspinup + else + spinup = 1 + endif + + + do it = 1,spinup + +!Step 1 - Initialize variables to 0 and extract the halo + + CA=0 + neighbours=0 + birth=0 + newlives=0 + neg=0 + newcell=0 + oldlives=0 + newval=0 + frac=0 + board_halo=0 + Wpert_halo=0 + Cpert_halo=0 + field_in=0 + maxlives = 0 + maxliveshigh =0 + livesout = 0 + + +!Step 4 - Compute the neighbourhood + + !Count the number of neighbours where perturbed massflux is larger than + !a threshold + + + ! do j=1,nyc + ! do i=1,nxc + ! ilives(i,j,nf)=int(real(nlives)*1.5*NOISE_B(i,j)) + ! if(vertvelhigh(i,j) > 0.)then ! .or. condition(i,j)<=0.)then + ! ilives(i,j,nf)=0 + ! endif + ! enddo + ! enddo + + do j=1,nyc + do i=1,nxc + field_in(i+(j-1)*nxc,1)=board(i,j,nf) + enddo + enddo + + call atmosphere_scalar_field_halo(board_halo,halo,nxch,nych,k_in,field_in) + + + do j=1,nyc + do i=1,nxc + ih=i+halo + jh=j+halo + neighbours(i,j)=board_halo(ih-1,jh-1,1)+board_halo(ih-1,jh,1)+ & + board_halo(ih-1,jh+1,1)+board_halo(ih,jh+1,1)+board_halo(ih+1,jh+1,1)+& + board_halo(ih+1,jh,1)+board_halo(ih+1,jh-1,1)+board_halo(ih,jh-1,1) + enddo + enddo + +!endif +! Step 5 - Check rules; + + !birth + + do j=1,nyc + do i=1,nxc + ! if(Wpert(i,j) < thresh(i,j) .and. (neighbours(i,j) < 2 .or. neighbours(i,j) > 3))then + if(neighbours(i,j) == 3 .or. neighbours(i,j) ==2)then + birth(i,j)=1 + endif + enddo + enddo + + !death + do j=1,nyc + do i=1,nxc + if(neighbours(i,j) < 2 .or. neighbours(i,j) > 3)then + !if(Wpert(i,j) > thresh(i,j) .and. (neighbours(i,j) < 2 .or. neighbours(i,j) > 3))then + lives(i,j,nf)=lives(i,j,nf) - 1 + endif + enddo + enddo + + do j=1,nyc + do i=1,nxc + if(lives(i,j,nf) < 0)then + lives(i,j,nf)=0 + endif + enddo + enddo + + do j=1,nyc + do i=1,nxc + if(birth(i,j)==1 .and. lives(i,j,nf)==0)then + newcell(i,j)=1 + endif + enddo + enddo + + do j=1,nyc + do i=1,nxc + lives(i,j,nf)=lives(i,j,nf)+newcell(i,j)*ilives(i,j,nf) !int(real(nlives)*1.5*NOISE_B(i,j)) + enddo + enddo + + do j=1,nyc + do i=1,nxc + if(neighbours(i,j)==3 .or. (board(i,j,nf)==1 .and. neighbours(i,j)==2))then + board(i,j,nf)=1 + else + board(i,j,nf)=0 + endif + enddo + enddo + + +enddo !spinup + +!COARSE-GRAIN BACK TO NWP GRID + + inci=ncells + incj=ncells + sub=ncells-1 + DO j=1,nlat + DO i=1,nlon + CA(i,j)=(SUM(lives(inci-sub:inci,incj-sub:incj,nf)))/real(ncells*ncells) + inci=inci+ncells + ENDDO + inci=ncells + incj=incj+ncells + ENDDO + +if(nca_plumes == .true.) 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 +!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 + endif + enddo + enddo + + kend=ceiling((ncells*ncells)/2.) + if (.not. allocated(V))then + allocate(V(kend)) + endif + if (.not. allocated(L))then + allocate(L(kend)) + endif + + ca_plumes(:,:)=0 + inci=ncells + incj=ncells + sub=ncells-1 + DO j=1,nlat + DO i=1,nlon + onegrid(1:ncells,1:ncells)=temp(inci-sub:inci,incj-sub:incj) + call plumes(V,L,onegrid,ncells,ncells,kend) + do k=1,kend + if (V(k) == 1)then + ca_plumes(i,j)=ca_plumes(i,j)+L(k) + endif + enddo + inci=inci+ncells + ENDDO + inci=ncells + incj=incj+ncells + ENDDO + +else + +ca_plumes(:,:)=0. + +endif ! nca_plumes + +end subroutine update_cells_sgs + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini_g,ilives_g, & + nlives,ncells,nfracseed,nseed,nthresh,ca_global, & + ca_sgs,nspinup,condition,vertvelhigh,nf,nca_plumes) + +implicit none + +integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca +integer, intent(in) :: iini_g(nxc,nyc,nca), ilives_g(nxc,nyc) +real, intent(in) :: condition(nxc,nyc),vertvelhigh(nxc,nyc) +real, intent(out) :: CA(nlon,nlat) +integer,intent(out) :: ca_plumes(nlon,nlat) +integer, intent(in) :: nlives, ncells, nseed, nspinup, nf +real, intent(in) :: nfracseed, nthresh +logical,intent(in) :: nca_plumes +real, dimension(nlon,nlat) :: frac +integer,allocatable,save :: board_g(:,:,:), lives_g(:,:,:) +integer,allocatable :: V(:),L(:) +integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize +integer :: haloh, ih, jh,lives_max,kend +real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh +real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) +integer, dimension(nxc,nyc) :: neighbours, birth, newlives, thresh +integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp,newseed +integer, dimension(ncells,ncells) :: onegrid +real,dimension(nxc,nyc) :: normlives +logical :: ca_global, ca_sgs +real :: Wpert(nxc,nyc),Cpert(nxc,nyc) + +!SGS parameters: +integer(8) :: count, count_rate, count_max, count_trunc +integer(8) :: iscale = 10000000000 +integer :: count5, count6 +type(random_stat) :: rstate +real :: dt, timescale, sigma, skew, kurt, acorr, gamma +real :: E_sq2, E2, g2, B_sq2, B2, sqrtdt,flamx2, tmp1, tmp1a +real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D +real, dimension(nxc*nyc) :: noise1D2, noise1D1 +real, allocatable, save :: sgs1(:,:,:),sgs2(:,:,:) +integer, dimension(nxch,nych,1) :: M_halo + + +!------------------------------------------------------------------------------------------------- +halo=1 +isize=nlon+2*halo +jsize=nlat+2*halo +k_in=1 + +ca_plumes=0. + + if (.not. allocated(board_g))then + allocate(board_g(nxc,nyc,nca)) + endif + if (.not. allocated(lives_g))then + allocate(lives_g(nxc,nyc,nca)) + endif + if (.not. allocated(sgs1))then + allocate(sgs1(nxc,nyc,nca)) + endif + if (.not. allocated(sgs2))then + allocate(sgs2(nxc,nyc,nca)) + endif + 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)) + endif + if(.not. allocated(Wpert_halo))then + allocate(Wpert_halo(nxch,nych,1)) + endif + if(.not. allocated(Cpert_halo))then + allocate(Cpert_halo(nxch,nych,1)) + endif + + !random numbers: + 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 + + call random_setseed(count5) + call random_number(noise1D1) + + call random_setseed(count6) + call random_number(noise1D2) + + !Put on 2D: + do j=1,nyc + do i=1,nxc + NOISE_A(i,j)=noise1D1(i+(j-1)*nxc)-0.5 + NOISE_B(i,j)=noise1D2(i+(j-1)*nxc) + enddo + enddo + + + if(kstep == 0)then + + do j=1,nyc + do i=1,nxc + board_g(i,j,nf) = iini_g(i,j,nf) + lives_g(i,j,nf) = ilives_g(i,j)*iini_g(i,j,nf) + enddo + enddo + + endif + + !Seed with new CA cells at each nseed step + newseed=0 + + if(mod(kstep,nseed) == 0 .and. kstep > 0.)then + do j=1,nyc + do i=1,nxc + if(board_g(i,j,nf) == 0 .and. NOISE_B(i,j)>0.95 )then + newseed(i,j)=1 + endif + board_g(i,j,nf) = board_g(i,j,nf) + newseed(i,j) + enddo + enddo + endif + + if(kstep == 0)then + spinup=nspinup + else + spinup = 1 + endif + + +do it=1,spinup + + +!Step 2 - Initialize variables to 0 and extract the halo + + neighbours=0 + birth=0 + newlives=0 + neg=0 + newcell=0 + oldlives=0 + newval=0 + frac=0 + board_halo=0 + Wpert_halo=0 + Cpert_halo=0 + field_in=0 + +!The input to scalar_field_halo needs to be 1D. +!take the updated board_g fields and extract the halo +! in order to have updated values in the halo region. + + + do j=1,nyc + do i=1,nxc + field_in(i+(j-1)*nxc,1)=board_g(i,j,nf) + enddo + enddo + + call atmosphere_scalar_field_halo(board_halo,halo,nxch,nych,k_in,field_in) + + + do j=1,nyc + do i=1,nxc + ih=i+halo + jh=j+halo + neighbours(i,j)=board_halo(ih-1,jh-1,1)+board_halo(ih-1,jh,1)+ & + board_halo(ih-1,jh+1,1)+board_halo(ih,jh+1,1)+board_halo(ih+1,jh+1,1)+& + board_halo(ih+1,jh,1)+board_halo(ih+1,jh-1,1)+board_halo(ih,jh-1,1) + enddo + enddo + + + + do j=1,nyc + do i=1,nxc + if(neighbours(i,j)==2 .or. neighbours(i,j)==3)then + birth(i,j)=1 + endif + enddo + enddo + + + do j=1,nyc + do i=1,nxc + if(neighbours(i,j)<2 .or. neighbours(i,j)>3)then + lives_g(i,j,nf)=lives_g(i,j,nf) - 1 + endif + enddo + enddo + + do j=1,nyc + do i=1,nxc + if(lives_g(i,j,nf)<0)then + lives_g(i,j,nf)=0 + endif + enddo + enddo + + do j=1,nyc + do i=1,nxc + if(birth(i,j)==1 .and. lives_g(i,j,nf)==0)then + newcell(i,j)=1 + endif + enddo + enddo + + + do j=1,nyc + do i=1,nxc + lives_g(i,j,nf)=lives_g(i,j,nf)+newcell(i,j)*ilives_g(i,j) + enddo + enddo + + + do j=1,nyc + do i=1,nxc + if( (board_g(i,j,nf) ==1 .and. (neighbours(i,j)==3 .or. neighbours(i,j)==2) ).or. (board_g(i,j,nf)==0 .and. neighbours(i,j)==3) )then + board_g(i,j,nf)=1 + else + board_g(i,j,nf)=0 + endif + enddo + enddo + +enddo !spinup + +!COARSE-GRAIN BACK TO NWP GRID + + inci=ncells + incj=ncells + sub=ncells-1 + DO j=1,nlat + DO i=1,nlon + CA(i,j)=(SUM(lives_g(inci-sub:inci,incj-sub:incj,nf)))/(ncells*ncells) + inci=inci+ncells + ENDDO + inci=ncells + incj=incj+ncells + ENDDO + +!lives_max=maxval(ilives_g) +!call mp_reduce_max(lives_max) +!CA(:,:) = (frac(:,:)/real(nlives/2)) + + +end subroutine update_cells_global + +end module update_ca diff --git a/update_ca.f90 b/update_ca.f90 deleted file mode 100644 index 10014d3..0000000 --- a/update_ca.f90 +++ /dev/null @@ -1,586 +0,0 @@ -module update_ca - -use atmosphere_mod, only: atmosphere_scalar_field_halo -use mersenne_twister, only: random_setseed,random_gauss,random_stat -use fv_mp_mod, only : mp_reduce_sum,mp_bcst,mp_reduce_min,mp_reduce_max - -implicit none - -!L. Bengtsson 2017-06 -!Evolve the cellular automata in time - - -contains - -subroutine 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) - -implicit none - -integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca -integer, intent(in) :: iini(nxc,nyc,nca), ilives(nxc,nyc) -real, intent(in) :: condition(nxc,nyc),vertvelhigh(nxc,nyc) -real, intent(out) :: CA(nlon,nlat) -integer,intent(out) :: ca_plumes(nlon,nlat) -integer, intent(in) :: nlives, ncells, nseed, nspinup, nf -real, intent(in) :: nfracseed, nthresh -logical,intent(in) :: nca_plumes -real, dimension(nlon,nlat) :: frac -integer,allocatable,save :: board(:,:,:), lives(:,:,:) -integer,allocatable :: V(:),L(:) -integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize -integer :: haloh, ih, jh,lives_max,kend -real :: thresh,threshc,threshk,wp_max,wp_min,mthresh,kthresh -real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) -integer, dimension(nxc,nyc) :: neighbours, birth, newlives -integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp -integer, dimension(ncells,ncells) :: onegrid -real,dimension(nxc,nyc) :: normlives -logical :: ca_global, ca_sgs -real :: Wpert(nxc,nyc),Cpert(nxc,nyc) - -!SGS parameters: -integer(8) :: count, count_rate, count_max, count_trunc -integer(8) :: iscale = 10000000000 -integer :: count5, count6 -type(random_stat) :: rstate -real :: dt, timescale, sigma, skew, kurt, acorr, gamma -real :: E_sq2, E2, g2, B_sq2, B2, sqrtdt,flamx2, tmp1, tmp1a -real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D -real, dimension(nxc*nyc) :: noise1D2, noise1D1 -real, allocatable, save :: sgs1(:,:,:),sgs2(:,:,:) -integer, dimension(nxch,nych,1) :: M_halo - - -!------------------------------------------------------------------------------------------------- -halo=1 -isize=nlon+2*halo -jsize=nlat+2*halo -k_in=1 - - if (.not. allocated(board))then - allocate(board(nxc,nyc,nca)) - endif - if (.not. allocated(lives))then - allocate(lives(nxc,nyc,nca)) - endif - if (.not. allocated(sgs1))then - allocate(sgs1(nxc,nyc,nca)) - endif - if (.not. allocated(sgs2))then - allocate(sgs2(nxc,nyc,nca)) - endif - 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)) - endif - if(.not. allocated(Wpert_halo))then - allocate(Wpert_halo(nxch,nych,1)) - endif - if(.not. allocated(Cpert_halo))then - allocate(Cpert_halo(nxch,nych,1)) - endif - - 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) - enddo - enddo - - endif - - if(kstep == 2)then - spinup=nspinup - else - spinup = 1 - endif - - else !ca_global - - if(kstep == 0)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) - enddo - enddo - - endif - - !Seed with new CA cells at each nseed step - if(mod(kstep,nseed) == 0)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) - enddo - enddo - - endif - - if(kstep == 0)then - spinup=nspinup - 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) then -!Compute the SGS and perturb the vertical velocity field: -!Read these values in from namelist, guided from LES data: - -dt=900. -timescale=21600.0 -sigma=0.8 -skew=0.8 -kurt=2.0 -acorr=exp(-dt/timescale) -gamma=-log(acorr)/dt - -!calculate coeffients for SGS with auto-correlation and use -!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) -if(kstep>=2)then -g2D=g2D+E2*condition(:,:) -endif -B_sq2=2.0*sigma**2*(gamma-E_sq2/2.0-(gamma-E_sq2)**2*skew**2/(8.0*E_sq2)) -B2=sqrt(B_sq2) -B2=0. - -sqrtdt=sqrt(dt) -flamx2=0.5*E_sq2+gamma - -endif - -do it=1,spinup - -if (ca_sgs) then - - !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 - call mp_bcst(count5) - call mp_bcst(count6) - - call random_setseed(count5,rstate) - call random_gauss(noise1D1,rstate) - - call random_setseed(count6,rstate) - call random_gauss(noise1D2,rstate) - !Put on 2D: - do j=1,nyc - do i=1,nxc - NOISE_A(i,j)=noise1D1(i+(j-1)*nxc) - NOISE_B(i,j)=noise1D2(i+(j-1)*nxc) - enddo - enddo - -tmp1=0. -tmp1a=0. -Wpert=0. -Cpert=0. - -if(kstep == 0)then - do j=1,nyc - do i=1,nxc - sgs1(i,j,nf)=sigma*NOISE_A(i,j) - sgs2(i,j,nf)=sigma*NOISE_B(i,j) - Cpert(i,j)=sgs1(i,j,nf) - Wpert(i,j)=sgs2(i,j,nf) - enddo - enddo - -else - - do j=1,nyc - do i=1,nxc - tmp1=sgs1(i,j,nf)-(flamx2*sgs1(i,j,nf)+0.5*E2*g2D(i,j))*dt + (B2*NOISE_A(i,j)*sqrtdt) & - + (g2D(i,j) + E2 * sgs1(i,j,nf))* NOISE_B(i,j)*sqrtdt - tmp1a=(tmp1+sgs1(i,j,nf))*0.5 - sgs1(i,j,nf)=sgs1(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 - - Cpert(i,j)=condition(i,j)*(1.0 + sgs1(i,j,nf)) - - enddo - enddo - - do j=1,nyc - do i=1,nxc - tmp1=sgs2(i,j,nf)-(flamx2*sgs2(i,j,nf)+0.5*E2*g2D(i,j))*dt + (B2*NOISE_A(i,j)*sqrtdt) & - + (g2D(i,j) + E2 * sgs2(i,j,nf))* NOISE_B(i,j)*sqrtdt - tmp1a=(tmp1+sgs2(i,j,nf))*0.5 - - 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)) - enddo - enddo - -endif - -endif !ca sgs true - - -!Step 2 - Initialize variables to 0 and extract the halo - - neighbours=0 - birth=0 - newlives=0 - neg=0 - newcell=0 - oldlives=0 - newval=0 - frac=0 - board_halo=0 - Wpert_halo=0 - Cpert_halo=0 - field_in=0 - -!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)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)then - field_in=0 - do j=1,nyc - do i=1,nxc - field_in(i+(j-1)*nxc,1)=Wpert(i,j) - enddo - enddo - - call atmosphere_scalar_field_halo(Wpert_halo,halo,nxch,nych,k_in,field_in) - - field_in=0 - do j=1,nyc - do i=1,nxc - field_in(i+(j-1)*nxc,1)=Cpert(i,j) - enddo - enddo - - call atmosphere_scalar_field_halo(Cpert_halo,halo,nxch,nych,k_in,field_in) - endif !sgs - - - - -!Step 4 - Compute the neighbourhood -if(ca_sgs)then !SGSmethod - - !Count the number of neighbours where perturbed massflux is larger than - !a threshold - - -if(nf==1)then !Deep convection - M_halo = 0 - do j=1,nych - do i=1,nxch - if(Wpert_halo(i,j,1) < nthresh .and. Cpert_halo(i,j,1) > nthresh)then - M_halo(i,j,1) = 1 - endif - enddo - enddo - elseif(nf==2)then !Shallow convection - M_halo=0 - do j=1,nych - do i=1,nxch - if(Wpert_halo(i,j,1) < nthresh .and. Cpert_halo(i,j,1)>nthresh)then - M_halo(i,j,1) = 1 - endif - enddo - enddo - elseif(nf==3)then !Turbulence - M_halo = 0 - do j=1,nych - do i=1,nxch - if(Wpert_halo(i,j,1) < nthresh .and. Cpert_halo(i,j,1) > nthresh)then - M_halo(i,j,1) = 1 - endif - enddo - enddo - elseif(nf==4)then !Radiation - M_halo = 0 - do j=1,nych - do i=1,nxch - if(Cpert_halo(i,j,1)>nthresh)then - M_halo(i,j,1) = 1 - endif - enddo - enddo - else !nf=5 Microphysics - M_halo = 0 - do j=1,nych - do i=1,nxch - if(Wpert_halo(i,j,1) < nthresh .and. Cpert_halo(i,j,1)>nthresh)then - M_halo(i,j,1) = 1 - endif - enddo - enddo -endif !nf - - do j=1,nyc - do i=1,nxc - ih=i+halo - jh=j+halo - neighbours(i,j)=M_halo(ih-1,jh-1,1)+M_halo(ih-1,jh,1)+ & - M_halo(ih-1,jh+1,1)+M_halo(ih,jh+1,1)+M_halo(ih+1,jh+1,1)+& - M_halo(ih+1,jh,1)+M_halo(ih+1,jh-1,1)+M_halo(ih,jh-1,1) - enddo - enddo - - -!CA stand alone method -else !global - - do j=1,nyc - do i=1,nxc - ih=i+halo - jh=j+halo - neighbours(i,j)=board_halo(ih-1,jh-1,1)+board_halo(ih-1,jh,1)+ & - board_halo(ih-1,jh+1,1)+board_halo(ih,jh+1,1)+board_halo(ih+1,jh+1,1)+& - board_halo(ih+1,jh,1)+board_halo(ih+1,jh-1,1)+board_halo(ih,jh-1,1) - enddo - enddo - -endif !sgs/global -! Step 5 - Check rules; the birth condition differs between SGS and GOL method - -if(ca_sgs)then !SGS - - if(nf==1)then - do j=1,nyc - do i=1,nxc - if((Wpert(i,j) < nthresh .and. Cpert_halo(i,j,1) > nthresh) .or. neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - elseif(nf==2)then - do j=1,nyc - do i=1,nxc - if((Wpert(i,j) < nthresh .and. Cpert_halo(i,j,1) > nthresh) .or. neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - elseif(nf==3)then - do j=1,nyc - do i=1,nxc - if((Wpert(i,j) < nthresh .and. Cpert_halo(i,j,1) > nthresh) .or. neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - elseif(nf==4)then - do j=1,nyc - do i=1,nxc - if(Cpert(i,j) > nthresh .or. neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - else !nf=5 - do j=1,nyc - do i=1,nxc - if((Wpert(i,j) < nthresh .and. Cpert_halo(i,j,1) > nthresh) .or. neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - endif - -else !GOL - - do j=1,nyc - do i=1,nxc - if(neighbours(i,j)==2 .or. neighbours(i,j)==3)then - birth(i,j)=1 - endif - enddo - enddo - -endif - - do j=1,nyc - do i=1,nxc - if(neighbours(i,j).ne.4 .or. neighbours(i,j).ne.5)then - lives(i,j,nf)=lives(i,j,nf) - 1 - endif - enddo - enddo - - do j=1,nyc - do i=1,nxc - if(lives(i,j,nf)<0)then - lives(i,j,nf)=0 - endif - enddo - enddo - - do j=1,nyc - do i=1,nxc - if(birth(i,j)==1 .and. lives(i,j,nf)==0)then - newcell(i,j)=1 - endif - enddo - enddo - - 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) - enddo - enddo - else !GOL - do j=1,nyc - do i=1,nxc - lives(i,j,nf)=lives(i,j,nf)+newcell(i,j)*ilives(i,j)!*nlives - enddo - enddo - endif - - do j=1,nyc - do i=1,nxc - if(neighbours(i,j)==3 .or. (board(i,j,nf)==1 .and. neighbours(i,j)==2))then - board(i,j,nf)=1 - else - board(i,j,nf)=0 - endif - enddo - enddo - -enddo !spinup - -!COARSE-GRAIN BACK TO NWP GRID - - inci=ncells - incj=ncells - sub=ncells-1 - DO j=1,nlat - DO i=1,nlon - frac(i,j)=(SUM(lives(inci-sub:inci,incj-sub:incj,nf)))/(ncells*ncells) - inci=inci+ncells - ENDDO - inci=ncells - incj=incj+ncells - ENDDO - -lives_max=maxval(ilives) -call mp_reduce_max(lives_max) - - if(ca_sgs)then - CA(:,:) = (frac(:,:)/lives_max) - else !global - CA(:,:) = (frac(:,:)/real(nlives)) - endif - - -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 -!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 - endif - enddo - enddo - - kend=ceiling((ncells*ncells)/2.) - if (.not. allocated(V))then - allocate(V(kend)) - endif - if (.not. allocated(L))then - allocate(L(kend)) - endif - - ca_plumes(:,:)=0 - inci=ncells - incj=ncells - sub=ncells-1 - DO j=1,nlat - DO i=1,nlon - onegrid(1:ncells,1:ncells)=temp(inci-sub:inci,incj-sub:incj) - call plumes(V,L,onegrid,ncells,ncells,kend) - do k=1,kend - if (V(k) == 1)then - ca_plumes(i,j)=ca_plumes(i,j)+L(k) - endif - enddo - inci=inci+ncells - ENDDO - inci=ncells - incj=incj+ncells - ENDDO - -else - -ca_plumes(:,:)=0. - -endif ! nca_plumes - -end subroutine update_cells - -end module update_ca From 41abe85d32a55013bca92e3425f7f0248b6e14d9 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Thu, 23 Jan 2020 21:22:19 +0000 Subject: [PATCH 02/10] ca updates global pattern --- cellular_automata_global.F90 | 124 ++++++++++++----------------------- 1 file changed, 42 insertions(+), 82 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index 1f67a76..8c47020 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -1,6 +1,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & nca,ncells,nlives,nfracseed,nseed,nthresh,ca_global,ca_sgs,iseed_ca, & - ca_smooth,nsmooth,nspinup,blocksize) + ca_smooth,nspinup,blocksize,nsmooth,ca_amplitude) use machine use update_ca, only: update_cells_sgs, update_cells_global @@ -27,7 +27,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !the flag ca_global integer,intent(in) :: kstep,ncells,nca,nlives,nseed,iseed_ca,nspinup,nsmooth -real,intent(in) :: nfracseed,nthresh +real,intent(in) :: nfracseed,nthresh,ca_amplitude logical,intent(in) :: ca_global, ca_sgs, ca_smooth integer, intent(in) :: nblks,nlev,blocksize type(GFS_coupling_type),intent(inout) :: Coupling(nblks) @@ -76,8 +76,8 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & ! WRITE(*,*)'Entering cellular automata calculations' ! Some security checks for namelist combinations: - if(nca > 5)then - write(0,*)'Namelist option nca cannot be larger than 5 - exiting' + if(nca > 3)then + write(0,*)'Namelist option nca cannot be larger than 3 - exiting' stop endif @@ -275,8 +275,17 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & ca_sgs,nspinup, condition, vertvelhigh,nf,nca_plumes) + + + + CA2 = CA2 + CA + +enddo !nca + +CA1 = CA if (ca_smooth) then + do nn=1,nsmooth !number of itterations for the smoothing. field_in=0. @@ -284,7 +293,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !get halo do j=1,nlat do i=1,nlon - field_in(i+(j-1)*nlon,1)=CA(i,j) + field_in(i+(j-1)*nlon,1)=CA1(i,j) enddo enddo @@ -296,33 +305,31 @@ subroutine cellular_automata_global(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)=(2.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)+& - 2.0*field_out(ih+1,jh-1,1))/20. + 2.0*field_out(ih+1,jh-1,1))/18. enddo enddo do j=1,nlat do i=1,nlon - CA(i,j)=field_smooth(i,j) + CA1(i,j)=field_smooth(i,j) enddo enddo enddo !nn -endif - - if(nf==1)then - CA1(:,:)=CA(:,:) - elseif(nf==2)then - CA2(:,:)=CA(:,:) - else - CA3(:,:)=CA(:,:) - endif +endif !smooth +! if(nf==1)then +! CA1(:,:)=CA(:,:) +! elseif(nf==2)then +! CA2(:,:)=CA(:,:) +! else +! CA3(:,:)=CA(:,:) +! endif -enddo !nca loop !!!!Post processing, should be made into a separate subroutine @@ -330,6 +337,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !!!!Construct linear combinations of CA1, CA2 and CA3 +!if (nf==1) then !Use min-max method to normalize range !!!CA1 Detmax(1)=maxval(CA1) @@ -343,6 +351,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo +!mean: CAmean=0. psum=0. csum=0. @@ -358,89 +367,40 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & CAmean=psum/csum -do j=1,nlat - do i=1,nlon - CA1(i,j)=(CA1(i,j)-CAmean)+1.0 !Can we compute the global median? - enddo -enddo - - - -!!!CA2 -Detmax(2)=maxval(CA2) -call mp_reduce_max(Detmax(2)) -Detmin(2)=minval(CA2) -call mp_reduce_min(Detmin(2)) - -do j=1,nlat - do i=1,nlon - CA2(i,j) = ((CA2(i,j) - Detmin(2))/(Detmax(2)-Detmin(2))) - enddo -enddo - -CAmean=0. -psum=0. -csum=0. -do j=1,nlat - do i=1,nlon - psum=psum+(CA2(i,j)) - csum=csum+1 - enddo -enddo - -call mp_reduce_sum(psum) -call mp_reduce_sum(csum) - -CAmean=psum/csum - -do j=1,nlat - do i=1,nlon - CA2(i,j)=(CA2(i,j)-CAmean)+1.0 !Can we compute the global median? - enddo -enddo +!std: +CAstdv=0. +sq_diff = 0. +do j=1,nlat + do i=1,nlon + sq_diff = sq_diff + (CA1(i,j)-CAmean)**2.0 + enddo +enddo +call mp_reduce_sum(sq_diff) -!!!CA3 -Detmax(3)=maxval(CA3) -call mp_reduce_max(Detmax(3)) -Detmin(3)=minval(CA3) -call mp_reduce_min(Detmin(3)) +CAstdv = sqrt(sq_diff/csum) -do j=1,nlat - do i=1,nlon - CA3(i,j) = ((CA3(i,j) - Detmin(3))/(Detmax(3)-Detmin(3))) - enddo -enddo +!Transform to mean of 1 and ca_amplitude standard deviation -CAmean=0. -psum=0. -csum=0. do j=1,nlat do i=1,nlon - psum=psum+(CA3(i,j)) - csum=csum+1 + CA1(i,j)=1.0 + (CA1(i,j)-CAmean)*(ca_amplitude/CAstdv) enddo enddo -call mp_reduce_sum(psum) -call mp_reduce_sum(csum) - -CAmean=psum/csum - do j=1,nlat do i=1,nlon - CA3(i,j)=(CA3(i,j)-CAmean)+1.0 !Can we compute the global median? + CA1(i,j)=min(max(CA1(i,j),0.),2.0) enddo enddo - !Put back into blocks 1D array to be passed to physics !or diagnostics output if(kstep < 1)then CA1(:,:)=1. -CA2(:,:)=1. -CA3(:,:)=1. +!CA2(:,:)=1. +!CA3(:,:)=1. endif From 72f81903f594139550df67b03910fde224f8b7b0 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 14 Feb 2020 16:36:03 +0000 Subject: [PATCH 03/10] latest CA updates SGS --- cellular_automata_global.F90 | 14 ++-- cellular_automata_sgs.F90 | 158 ++++++----------------------------- update_ca.F90 | 30 +++++-- 3 files changed, 54 insertions(+), 148 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index 8c47020..6f5bad0 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -46,7 +46,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:),Detfield(:,:,:) real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:) real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:),tmp(:,:) -real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),condition(:,:),rho(:,:),cape(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),condition(:,:),rho(:,:),conditiongrid(:,:) real(kind=kind_phys), allocatable :: noise1D(:),vertvelhigh(:,:),noise(:,:,:) real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,lambda real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt @@ -127,7 +127,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & allocate(ilives_g(nxc,nyc)) allocate(vertvelhigh(nxc,nyc)) allocate(condition(nxc,nyc)) - allocate(cape(nlon,nlat)) + allocate(conditiongrid(nlon,nlat)) allocate(Detfield(nlon,nlat,nca)) allocate(CA(nlon,nlat)) allocate(ca_plumes(nlon,nlat)) @@ -145,7 +145,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & cloud(:,:)=0. humidity(:,:)=0. condition(:,:)=0. - cape(:,:)=0. + conditiongrid(:,:)=0. vertvelhigh(:,:)=0. noise(:,:,:) = 0.0 noise1D(:) = 0.0 @@ -175,7 +175,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & 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) + conditiongrid(i,j) = Coupling(blk)%condition(ix) surfp(i,j) = Statein(blk)%pgr(ix) humidity(i,j)=Statein(blk)%qgrs(ix,13,1) !about 850 hpa do k = 1,29 !Lower troposphere: level 29 is about 350hPa @@ -276,10 +276,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & - - - CA2 = CA2 + CA - enddo !nca CA1 = CA @@ -421,7 +417,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & deallocate(pressure) deallocate(humidity) deallocate(dp) - deallocate(cape) + deallocate(conditiongrid) deallocate(rho) deallocate(surfp) deallocate(vertvelmean) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index b0a7822..56fc169 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -57,11 +57,11 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:),Detfield(:,:,:) real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:),uwind(:,:) real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:),shalp(:,:),gamt(:,:) -real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),rho(:,:),cape(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),rho(:,:),conditiongrid(:,:) real(kind=kind_phys), allocatable :: CA_DEEP(:,:),CA_TURB(:,:),CA_SHAL(:,:) real(kind=kind_phys), allocatable :: noise1D(:),vertvelhigh(:,:),noise(:,:,:) real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,lambda -real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt +real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt,condmax logical,save :: block_message=.true. logical :: nca_plumes @@ -80,7 +80,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & halo=1 k_in=1 -nca_plumes = .false. +nca_plumes = .true. !---------------------------------------------------------------------------- ! Get information about the compute domain, allocate fields on this ! domain @@ -105,13 +105,14 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & nych=nyc+2*halo !Allocate fields: + levs=nlev allocate(cloud(nlon,nlat)) - allocate(omega(nlon,nlat,29)) - allocate(pressure(nlon,nlat,29)) + allocate(omega(nlon,nlat,levs)) + allocate(pressure(nlon,nlat,levs)) allocate(humidity(nlon,nlat)) allocate(uwind(nlon,nlat)) - allocate(dp(nlon,nlat,29)) + allocate(dp(nlon,nlat,levs)) allocate(rho(nlon,nlat)) allocate(surfp(nlon,nlat)) allocate(vertvelmean(nlon,nlat)) @@ -125,7 +126,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & allocate(ilives_g(nxc,nyc)) allocate(vertvelhigh(nxc,nyc)) allocate(condition(nxc,nyc)) - allocate(cape(nlon,nlat)) + allocate(conditiongrid(nlon,nlat)) allocate(shalp(nlon,nlat)) allocate(gamt(nlon,nlat)) allocate(Detfield(nlon,nlat,nca)) @@ -145,7 +146,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & humidity(:,:)=0. uwind(:,:) = 0. condition(:,:)=0. - cape(:,:)=0. + conditiongrid(:,:)=0. vertvelhigh(:,:)=0. ca_plumes(:,:) = 0 noise(:,:,:) = 0.0 @@ -175,7 +176,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & i = Atm_block%index(blk)%ii(ix) - isc + 1 j = Atm_block%index(blk)%jj(ix) - jsc + 1 uwind(i,j) = Statein(blk)%ugrs(ix,30) - cape(i,j) = Coupling(blk)%cape(ix) + conditiongrid(i,j) = Coupling(blk)%condition(ix) shalp(i,j) = Coupling(blk)%ca_shal(ix) gamt(i,j) = Coupling(blk)%ca_turb(ix) surfp(i,j) = Statein(blk)%pgr(ix) @@ -267,7 +268,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) !cape(inci/ncells,incj/ncells) + condition(i,j)=conditiongrid(inci/ncells,incj/ncells) !conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -278,9 +279,12 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo + condmax=maxval(condition) + call mp_reduce_max(condmax) + do j = 1,nyc do i = 1,nxc - ilives(i,j,nf)=nlives + ilives(i,j,nf)=real(nlives)*(condition(i,j)/condmax) enddo enddo @@ -289,7 +293,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=shalp(inci/ncells,incj/ncells) + condition(i,j)=conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -300,40 +304,22 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo + condmax=maxval(condition) + call mp_reduce_max(condmax) + do j = 1,nyc do i = 1,nxc - ilives(i,j,nf)=nlives !int(real(nlives)*1.5*noise(i,j,nf)) + ilives(i,j,nf)=real(nlives)*(condition(i,j)/condmax) enddo enddo - elseif(nf==3)then - inci=ncells - incj=ncells - do j=1,nyc - do i=1,nxc - condition(i,j)=gamt(inci/ncells,incj/ncells) - if(i.eq.inci)then - inci=inci+ncells - endif - enddo - inci=ncells - if(j.eq.incj)then - incj=incj+ncells - endif - enddo + else - do j = 1,nyc - do i = 1,nxc - ilives(i,j,nf)=nlives - enddo - enddo - - elseif(nf==4)then inci=ncells incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=vertvelmean(inci/ncells,incj/ncells) + condition(i,j)=conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -344,35 +330,15 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo - do j = 1,nyc - do i = 1,nxc - ilives(i,j,nf)=int(real(nlives)*1.5*noise(i,j,nf)) - enddo - enddo - - else - inci=ncells - incj=ncells - do j=1,nyc - do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) - if(i.eq.inci)then - inci=inci+ncells - endif - enddo - inci=ncells - if(j.eq.incj)then - incj=incj+ncells - endif - enddo + condmax=maxval(condition) + call mp_reduce_max(condmax) do j = 1,nyc do i = 1,nxc - ilives(i,j,nf)=int(real(nlives)*1.5*noise(i,j,nf)) + ilives(i,j,nf)=real(nlives)*(condition(i,j)/condmax) enddo enddo - endif !nf @@ -418,16 +384,6 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & if(kstep > 1)then -!Box-Cox method to make output range less skewed -lambda=0.38 -do j=1,nlat - do i=1,nlon - if(CA_DEEP(i,j).NE.0.)then - CA_DEEP(i,j)=((CA_DEEP(i,j)**lambda)-1.0)/lambda - endif - enddo -enddo - !Use min-max method to normalize range Detmax(1)=maxval(CA_DEEP,CA_DEEP.NE.0.) call mp_reduce_max(Detmax(1)) @@ -472,28 +428,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & Detmin(1) = minval(CA_DEEP,CA_DEEP.NE.0) call mp_reduce_min(Detmin(1)) -!Set range outside positive vertical mean vertical velocity to 0. -!do j=1,nlat -! do i=1,nlon -! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward -! CA_DEEP(i,j)=0. -! endif -! enddo -!enddo - - - !Shallow convection ============================================================ -!Box-Cox method to make range less skewed - -lambda=0.38 -do j=1,nlat - do i=1,nlon - if(CA_SHAL(i,j).NE.0.)then - CA_SHAL(i,j)=((CA_SHAL(i,j)**lambda)-1.0)/lambda - endif - enddo -enddo !Use min-max method to normalize range Detmax(2)=maxval(CA_SHAL,CA_SHAL.NE.0) @@ -535,27 +470,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo -!Set range outside positive vertical mean vertical velocity to 0. -!do j=1,nlat -! do i=1,nlon -! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward ! -! CA_SHAL(i,j)=0. -! endif -! enddo -!enddo - - !Turbulence ============================================================================= -!Box-Cox method to make the range less skewed - -lambda=0.38 -do j=1,nlat - do i=1,nlon - if(CA_TURB(i,j).NE.0.)then - CA_TURB(i,j)=((CA_TURB(i,j)**lambda)-1.0)/lambda - endif - enddo -enddo !Use min-max method to normalize range Detmax(3)=maxval(CA_TURB,CA_TURB.NE.0) @@ -597,27 +512,8 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo -!Set range outside positive vertical mean vertical velocity to 0. -!do j=1,nlat -! do i=1,nlon -! if(vertvelmean(i,j)>0.)then !larger than because units are Pa/s so negative wind is upward ! -! CA_TURB(i,j)=0. -! endif -! enddo -!enddo - endif !kstep >1 -!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 - ca_plumes(i,j)=20 - endif - enddo -enddo - !Put back into blocks 1D array to be passed to physics !or diagnostics output @@ -625,7 +521,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & 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 - Diag(blk)%ca_deep(ix)=CA_DEEP(i,j) + Diag(blk)%ca_deep(ix)=ca_plumes(i,j) Diag(blk)%ca_turb(ix)=CA_TURB(i,j) Diag(blk)%ca_shal(ix)=CA_SHAL(i,j) Coupling(blk)%ca_deep(ix)=CA_DEEP(i,j) @@ -639,7 +535,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & deallocate(pressure) deallocate(humidity) deallocate(dp) - deallocate(cape) + deallocate(conditiongrid) deallocate(shalp) deallocate(gamt) deallocate(rho) diff --git a/update_ca.F90 b/update_ca.F90 index 106434a..1bf9bfa 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -34,13 +34,14 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i real, dimension(nlon,nlat) :: frac integer, dimension(nlon,nlat) :: maxlives integer,allocatable,save :: board(:,:,:), lives(:,:,:) -integer,allocatable :: V(:),L(:) +integer,allocatable :: V(:),L(:),B(:) +integer,allocatable :: AG(:,:) integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize integer :: haloh, ih, jh,lives_max,kend real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) integer, dimension(nxc,nyc) :: neighbours, birth, newlives,thresh,maxliveshigh -integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp +integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp,newseed integer, dimension(ncells,ncells) :: onegrid real,dimension(nxc,nyc) :: livesout logical :: ca_global, ca_sgs @@ -159,11 +160,14 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i !Seed with new CA cells at each nseed step + newseed = 0 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,nf)*iini(i,j,nf) + if(board(i,j,nf) == 0 .and. NOISE_B(i,j)>0.95 )then + newseed(i,j) = 1 + endif + board(i,j,nf) = board(i,j,nf) + newseed(i,j) enddo enddo @@ -327,6 +331,12 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i if (.not. allocated(L))then allocate(L(kend)) endif + if (.not. allocated(B))then + allocate(B(kend)) + endif + if (.not. allocated(AG))then + allocate(AG(ncells,ncells)) + endif ca_plumes(:,:)=0 inci=ncells @@ -334,13 +344,17 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i sub=ncells-1 DO j=1,nlat DO i=1,nlon + B(:)=0 + L(:)=0 + V(:)=0 onegrid(1:ncells,1:ncells)=temp(inci-sub:inci,incj-sub:incj) - call plumes(V,L,onegrid,ncells,ncells,kend) + call plumes(V,L,AG,onegrid,ncells,ncells,kend) do k=1,kend - if (V(k) == 1)then - ca_plumes(i,j)=ca_plumes(i,j)+L(k) - endif + if(V(k)==1)then + B(k)=L(k) !to avoid considering clusters of 0 + endif enddo + ca_plumes(i,j)=MAXVAL(B(1:kend)) inci=inci+ncells ENDDO inci=ncells From 63768a79d0c2f1f4c0f9080d9bd9e651af90fcd9 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 18 Feb 2020 23:18:27 +0000 Subject: [PATCH 04/10] cleaning of final code --- cellular_automata_global.F90 | 161 ++++++++--------------------------- cellular_automata_sgs.F90 | 17 +++- plumes.f90 | 6 +- update_ca.F90 | 132 +++------------------------- 4 files changed, 63 insertions(+), 253 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index 6f5bad0..ff417f2 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -42,16 +42,14 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & integer :: blocksz,levs integer(8) :: count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 -integer, allocatable :: iini(:,:,:),ilives(:,:),iini_g(:,:,:),ilives_g(:,:),ca_plumes(:,:) -real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:),Detfield(:,:,:) -real(kind=kind_phys), allocatable :: omega(:,:,:),pressure(:,:,:),cloud(:,:),humidity(:,:) -real(kind=kind_phys), allocatable :: vertvelsum(:,:),vertvelmean(:,:),dp(:,:,:),surfp(:,:),tmp(:,:) -real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),condition(:,:),rho(:,:),conditiongrid(:,:) +integer, allocatable :: iini_g(:,:,:),ilives_g(:,:) +real(kind=kind_phys), allocatable :: field_out(:,:,:), field_in(:,:),field_smooth(:,:) +real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:) real(kind=kind_phys), allocatable :: noise1D(:),vertvelhigh(:,:),noise(:,:,:) -real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv,count1,lambda -real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca),phi,stdev,delt +real(kind=kind_phys) :: psum,csum,CAmean,sq_diff,CAstdv +real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca) logical,save :: block_message=.true. -logical :: nca_plumes + !nca :: switch for number of cellular automata to be used. !ca_global :: switch for global cellular automata @@ -63,12 +61,11 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & ! 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. + halo=1 k_in=1 -nca_plumes = .false. !---------------------------------------------------------------------------- ! Get information about the compute domain, allocate fields on this ! domain @@ -109,55 +106,30 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !Allocate fields: - allocate(cloud(nlon,nlat)) - allocate(omega(nlon,nlat,29)) - allocate(pressure(nlon,nlat,29)) - allocate(humidity(nlon,nlat)) - allocate(dp(nlon,nlat,29)) - allocate(rho(nlon,nlat)) - allocate(surfp(nlon,nlat)) - allocate(vertvelmean(nlon,nlat)) - allocate(vertvelsum(nlon,nlat)) allocate(field_in(nlon*nlat,1)) allocate(field_out(isize,jsize,1)) allocate(field_smooth(nlon,nlat)) - allocate(iini(nxc,nyc,nca)) - allocate(ilives(nxc,nyc)) allocate(iini_g(nxc,nyc,nca)) allocate(ilives_g(nxc,nyc)) allocate(vertvelhigh(nxc,nyc)) - allocate(condition(nxc,nyc)) - allocate(conditiongrid(nlon,nlat)) - allocate(Detfield(nlon,nlat,nca)) allocate(CA(nlon,nlat)) - allocate(ca_plumes(nlon,nlat)) allocate(CA1(nlon,nlat)) allocate(CA2(nlon,nlat)) allocate(CA3(nlon,nlat)) - allocate(tmp(nlon,nlat)) allocate(noise(nxc,nyc,nca)) allocate(noise1D(nxc*nyc)) !Initialize: - Detfield(:,:,:)=0. - vertvelmean(:,:) =0. - vertvelsum(:,:)=0. - cloud(:,:)=0. - humidity(:,:)=0. - condition(:,:)=0. - conditiongrid(:,:)=0. - vertvelhigh(:,:)=0. + noise(:,:,:) = 0.0 noise1D(:) = 0.0 - iini(:,:,:) = 0 - ilives(:,:) = 0 iini_g(:,:,:) = 0 ilives_g(:,:) = 0 CA1(:,:) = 0.0 CA2(:,:) = 0.0 CA3(:,:) = 0.0 - ca_plumes(:,:) = 0 + !Put the blocks of model fields into a 2d array levs=nlev blocksz=blocksize @@ -171,43 +143,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & jsc = Atm_block%jsc 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 - conditiongrid(i,j) = Coupling(blk)%condition(ix) - surfp(i,j) = Statein(blk)%pgr(ix) - humidity(i,j)=Statein(blk)%qgrs(ix,13,1) !about 850 hpa - do k = 1,29 !Lower troposphere: level 29 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 - enddo - enddo - -!Compute layer averaged vertical velocity (Pa/s) - vertvelsum=0. - vertvelmean=0. - do j=1,nlat - do i =1,nlon - dp(i,j,1)=(surfp(i,j)-pressure(i,j,1)) - do k=2,29 - dp(i,j,k)=(pressure(i,j,k-1)-pressure(i,j,k)) - enddo - count1=0. - do k=1,29 - count1=count1+1. - vertvelsum(i,j)=vertvelsum(i,j)+(omega(i,j,k)*dp(i,j,k)) - enddo - enddo - enddo - - do j=1,nlat - do i=1,nlon - vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,29)) - enddo - enddo - !Generate random number, following stochastic physics code: do nf=1,nca @@ -253,7 +188,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !we here set the "condition" variable to a different model field depending !on nf. (this is not used if ca_global = .true.) - do nf=1,nca !update each ca do j = 1,nyc @@ -270,16 +204,10 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & CA(:,:)=0. - call update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini_g,ilives_g, & - nlives, ncells, nfracseed, nseed,nthresh, ca_global, & - ca_sgs,nspinup, condition, vertvelhigh,nf,nca_plumes) - - + call update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,iini_g,ilives_g, & + nlives, ncells, nfracseed, nseed,nthresh, nspinup,nf) -enddo !nca - -CA1 = CA - + if (ca_smooth) then do nn=1,nsmooth !number of itterations for the smoothing. @@ -289,7 +217,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !get halo do j=1,nlat do i=1,nlon - field_in(i+(j-1)*nlon,1)=CA1(i,j) + field_in(i+(j-1)*nlon,1)=CA(i,j) enddo enddo @@ -311,39 +239,23 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon - CA1(i,j)=field_smooth(i,j) + CA(i,j)=field_smooth(i,j) enddo enddo enddo !nn endif !smooth -! if(nf==1)then -! CA1(:,:)=CA(:,:) -! elseif(nf==2)then -! CA2(:,:)=CA(:,:) -! else -! CA3(:,:)=CA(:,:) -! endif - - - !!!!Post processing, should be made into a separate subroutine -!!!!Construct linear combinations of CA1, CA2 and CA3 - - -!if (nf==1) then -!Use min-max method to normalize range -!!!CA1 -Detmax(1)=maxval(CA1) +Detmax(1)=maxval(CA) call mp_reduce_max(Detmax(1)) -Detmin(1)=minval(CA1) +Detmin(1)=minval(CA) call mp_reduce_min(Detmin(1)) do j=1,nlat do i=1,nlon - CA1(i,j) = ((CA1(i,j) - Detmin(1))/(Detmax(1)-Detmin(1))) + CA(i,j) = ((CA(i,j) - Detmin(1))/(Detmax(1)-Detmin(1))) enddo enddo @@ -353,7 +265,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & csum=0. do j=1,nlat do i=1,nlon - psum=psum+(CA1(i,j)) + psum=psum+(CA(i,j)) csum=csum+1 enddo enddo @@ -368,7 +280,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & sq_diff = 0. do j=1,nlat do i=1,nlon - sq_diff = sq_diff + (CA1(i,j)-CAmean)**2.0 + sq_diff = sq_diff + (CA(i,j)-CAmean)**2.0 enddo enddo @@ -380,25 +292,31 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon - CA1(i,j)=1.0 + (CA1(i,j)-CAmean)*(ca_amplitude/CAstdv) + CA(i,j)=1.0 + (CA(i,j)-CAmean)*(ca_amplitude/CAstdv) enddo enddo do j=1,nlat do i=1,nlon - CA1(i,j)=min(max(CA1(i,j),0.),2.0) + CA(i,j)=min(max(CA(i,j),0.),2.0) enddo enddo - !Put back into blocks 1D array to be passed to physics !or diagnostics output if(kstep < 1)then -CA1(:,:)=1. -!CA2(:,:)=1. -!CA3(:,:)=1. +CA(:,:)=1. endif + if(nf==1)then + CA1(:,:)=CA(:,:) + elseif(nf==2)then + CA2(:,:)=CA(:,:) + else + CA3(:,:)=CA(:,:) + endif + +enddo !nf do blk = 1, Atm_block%nblks do ix = 1,Atm_block%blksz(blk) @@ -413,24 +331,13 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo - deallocate(omega) - deallocate(pressure) - deallocate(humidity) - deallocate(dp) - deallocate(conditiongrid) - deallocate(rho) - deallocate(surfp) - deallocate(vertvelmean) - deallocate(vertvelsum) + deallocate(field_in) deallocate(field_out) deallocate(field_smooth) - deallocate(iini) - deallocate(ilives) - deallocate(condition) - deallocate(Detfield) + deallocate(iini_g) + deallocate(ilives_g) deallocate(CA) - deallocate(ca_plumes) deallocate(CA1) deallocate(CA2) deallocate(CA3) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index 56fc169..72b535f 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -364,8 +364,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & !If ca-global is used, then nca independent CAs are called and weighted together to create one field; CA call update_cells_sgs(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) + nlives, ncells, nfracseed, nseed,nthresh,nspinup,nf,nca_plumes) if(nf==1)then CA_DEEP(:,:)=CA(:,:) @@ -514,6 +513,16 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & endif !kstep >1 +do j=1,nlat + do i=1,nlon + if(conditiongrid(i,j) == 0)then + CA_DEEP(i,j)=0. + ca_plumes(i,j)=0. + endif + enddo +enddo + + !Put back into blocks 1D array to be passed to physics !or diagnostics output @@ -522,9 +531,9 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & i = Atm_block%index(blk)%ii(ix) - isc + 1 j = Atm_block%index(blk)%jj(ix) - jsc + 1 Diag(blk)%ca_deep(ix)=ca_plumes(i,j) - Diag(blk)%ca_turb(ix)=CA_TURB(i,j) + Diag(blk)%ca_turb(ix)=conditiongrid(i,j) Diag(blk)%ca_shal(ix)=CA_SHAL(i,j) - Coupling(blk)%ca_deep(ix)=CA_DEEP(i,j) + Coupling(blk)%ca_deep(ix)=ca_plumes(i,j) Coupling(blk)%ca_turb(ix)=CA_TURB(i,j) Coupling(blk)%ca_shal(ix)=CA_SHAL(i,j) enddo diff --git a/plumes.f90 b/plumes.f90 index 410cc47..4714870 100644 --- a/plumes.f90 +++ b/plumes.f90 @@ -1,4 +1,4 @@ -subroutine plumes(V,L,a,row,col,kend) +subroutine plumes(V,L,AG,a,row,col,kend) implicit none !!January 2018 adapted from Mathworkds "ISLANDS" @@ -22,10 +22,10 @@ subroutine plumes(V,L,a,row,col,kend) integer, intent(in) :: row,col,kend integer, intent(in), dimension(row,col) :: a -integer, intent(out), dimension(kend) :: V,L +integer, intent(out), dimension(kend) :: L,V integer :: cnt,pp,mm,kk integer :: i,j,cntr,gg,hh,ii,jj,idxx,IDX -integer, dimension(row,col) :: AG +integer, intent(out), dimension(row,col) :: AG diff --git a/update_ca.F90 b/update_ca.F90 index 1bf9bfa..77771ee 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -17,17 +17,15 @@ module update_ca contains subroutine update_cells_sgs(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) + nlives,ncells,nfracseed,nseed,nthresh,nspinup,nf,nca_plumes) implicit none integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca integer, intent(in) :: iini(nxc,nyc,nca) integer, intent(inout) :: ilives(nxc,nyc,nca) -real, intent(in) :: condition(nxc,nyc),vertvelhigh(nxc,nyc) real, intent(out) :: CA(nlon,nlat) -integer,intent(out) :: ca_plumes(nlon,nlat) +integer, intent(out) :: ca_plumes(nlon,nlat) integer, intent(in) :: nlives, ncells, nseed, nspinup, nf real, intent(in) :: nfracseed, nthresh logical,intent(in) :: nca_plumes @@ -37,27 +35,19 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i integer,allocatable :: V(:),L(:),B(:) integer,allocatable :: AG(:,:) integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize -integer :: haloh, ih, jh,lives_max,kend +integer :: ih, jh,kend real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh -real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) +real, allocatable :: field_in(:,:),board_halo(:,:,:) integer, dimension(nxc,nyc) :: neighbours, birth, newlives,thresh,maxliveshigh integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp,newseed integer, dimension(ncells,ncells) :: onegrid -real,dimension(nxc,nyc) :: livesout -logical :: ca_global, ca_sgs -real :: Wpert(nxc,nyc),Cpert(nxc,nyc) -!SGS parameters: integer(8) :: count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 integer :: count5, count6 type(random_stat) :: rstate -real :: dt, timescale, sigma, skew, kurt, acorr, gamma -real :: E_sq2, E2, g2, B_sq2, B2, sqrtdt,flamx2, tmp1, tmp1a -real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D +real, dimension(nxc,nyc) :: NOISE_A, NOISE_B real, dimension(nxc*nyc) :: noise1D2, noise1D1 -real, allocatable, save :: sgs1(:,:,:),sgs2(:,:,:) -integer, dimension(nxch,nych,1) :: M_halo !------------------------------------------------------------------------------------------------- @@ -72,27 +62,13 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i if (.not. allocated(lives))then allocate(lives(nxc,nyc,nca)) endif - if (.not. allocated(sgs1))then - allocate(sgs1(nxc,nyc,nca)) - endif - if (.not. allocated(sgs2))then - allocate(sgs2(nxc,nyc,nca)) - endif 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)) endif - if(.not. allocated(Wpert_halo))then - allocate(Wpert_halo(nxch,nych,1)) - endif - if(.not. allocated(Cpert_halo))then - allocate(Cpert_halo(nxch,nych,1)) - endif - -! Some initial white noise generation: (could be SGS) -!Random seed for SGS + noise1D1 = 0.0 noise1D2 = 0.0 @@ -110,36 +86,6 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i call random_setseed(count6) call random_number(noise1D2) - !Put on 2D: -! do j=1,nyc -! do i=1,nxc -! NOISE_A(i,j)=noise1D1(i+(j-1)*nxc)-0.5 -! NOISE_B(i,j)=noise1D2(i+(j-1)*nxc) -! enddo -! enddo - -!Wpert=0. -!Cpert=0. - -! do j=1,nyc -! do i=1,nxc -! if(vertvelhigh(i,j) < 0. )then !upward motion -! Wpert(i,j)=vertvelhigh(i,j)*(1.0 + NOISE_A(i,j)) -! endif -! enddo -! enddo - -!endif - - -! do j=1,nyc -! do i=1,nxc -! if(vertvelhigh(i,j) < 0.)then !upward vertical motion -! thresh(i,j)=vertvelhigh(i,j) !condition(i,j) -! endif -! enddo -! enddo - if(kstep <= 1)then do j=1,nyc do i=1,nxc @@ -194,28 +140,12 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i newval=0 frac=0 board_halo=0 - Wpert_halo=0 - Cpert_halo=0 field_in=0 maxlives = 0 maxliveshigh =0 - livesout = 0 - + !Step 4 - Compute the neighbourhood - - !Count the number of neighbours where perturbed massflux is larger than - !a threshold - - - ! do j=1,nyc - ! do i=1,nxc - ! ilives(i,j,nf)=int(real(nlives)*1.5*NOISE_B(i,j)) - ! if(vertvelhigh(i,j) > 0.)then ! .or. condition(i,j)<=0.)then - ! ilives(i,j,nf)=0 - ! endif - ! enddo - ! enddo do j=1,nyc do i=1,nxc @@ -236,14 +166,12 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i enddo enddo -!endif ! Step 5 - Check rules; !birth do j=1,nyc do i=1,nxc - ! if(Wpert(i,j) < thresh(i,j) .and. (neighbours(i,j) < 2 .or. neighbours(i,j) > 3))then if(neighbours(i,j) == 3 .or. neighbours(i,j) ==2)then birth(i,j)=1 endif @@ -254,7 +182,6 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i do j=1,nyc do i=1,nxc if(neighbours(i,j) < 2 .or. neighbours(i,j) > 3)then - !if(Wpert(i,j) > thresh(i,j) .and. (neighbours(i,j) < 2 .or. neighbours(i,j) > 3))then lives(i,j,nf)=lives(i,j,nf) - 1 endif enddo @@ -278,7 +205,7 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i do j=1,nyc do i=1,nxc - lives(i,j,nf)=lives(i,j,nf)+newcell(i,j)*ilives(i,j,nf) !int(real(nlives)*1.5*NOISE_B(i,j)) + lives(i,j,nf)=lives(i,j,nf)+newcell(i,j)*ilives(i,j,nf) enddo enddo @@ -371,45 +298,32 @@ end subroutine update_cells_sgs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini_g,ilives_g, & - nlives,ncells,nfracseed,nseed,nthresh,ca_global, & - ca_sgs,nspinup,condition,vertvelhigh,nf,nca_plumes) +subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,iini_g,ilives_g, & + nlives,ncells,nfracseed,nseed,nthresh,nspinup,nf) implicit none integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca integer, intent(in) :: iini_g(nxc,nyc,nca), ilives_g(nxc,nyc) -real, intent(in) :: condition(nxc,nyc),vertvelhigh(nxc,nyc) real, intent(out) :: CA(nlon,nlat) -integer,intent(out) :: ca_plumes(nlon,nlat) integer, intent(in) :: nlives, ncells, nseed, nspinup, nf real, intent(in) :: nfracseed, nthresh -logical,intent(in) :: nca_plumes real, dimension(nlon,nlat) :: frac integer,allocatable,save :: board_g(:,:,:), lives_g(:,:,:) integer,allocatable :: V(:),L(:) integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize -integer :: haloh, ih, jh,lives_max,kend +integer :: ih, jh real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh -real, allocatable :: field_in(:,:),board_halo(:,:,:), Wpert_halo(:,:,:),Cpert_halo(:,:,:) +real, allocatable :: field_in(:,:),board_halo(:,:,:) integer, dimension(nxc,nyc) :: neighbours, birth, newlives, thresh integer, dimension(nxc,nyc) :: neg, newcell, oldlives, newval,temp,newseed -integer, dimension(ncells,ncells) :: onegrid -real,dimension(nxc,nyc) :: normlives -logical :: ca_global, ca_sgs -real :: Wpert(nxc,nyc),Cpert(nxc,nyc) -!SGS parameters: integer(8) :: count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 integer :: count5, count6 type(random_stat) :: rstate -real :: dt, timescale, sigma, skew, kurt, acorr, gamma -real :: E_sq2, E2, g2, B_sq2, B2, sqrtdt,flamx2, tmp1, tmp1a real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D real, dimension(nxc*nyc) :: noise1D2, noise1D1 -real, allocatable, save :: sgs1(:,:,:),sgs2(:,:,:) -integer, dimension(nxch,nych,1) :: M_halo !------------------------------------------------------------------------------------------------- @@ -418,32 +332,19 @@ subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plume jsize=nlat+2*halo k_in=1 -ca_plumes=0. - if (.not. allocated(board_g))then allocate(board_g(nxc,nyc,nca)) endif if (.not. allocated(lives_g))then allocate(lives_g(nxc,nyc,nca)) endif - if (.not. allocated(sgs1))then - allocate(sgs1(nxc,nyc,nca)) - endif - if (.not. allocated(sgs2))then - allocate(sgs2(nxc,nyc,nca)) - endif 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)) endif - if(.not. allocated(Wpert_halo))then - allocate(Wpert_halo(nxch,nych,1)) - endif - if(.not. allocated(Cpert_halo))then - allocate(Cpert_halo(nxch,nych,1)) - endif + !random numbers: noise1D1 = 0.0 @@ -515,8 +416,6 @@ subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plume newval=0 frac=0 board_halo=0 - Wpert_halo=0 - Cpert_halo=0 field_in=0 !The input to scalar_field_halo needs to be 1D. @@ -612,11 +511,6 @@ subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plume incj=incj+ncells ENDDO -!lives_max=maxval(ilives_g) -!call mp_reduce_max(lives_max) -!CA(:,:) = (frac(:,:)/real(nlives/2)) - - end subroutine update_cells_global end module update_ca From 42a6b47ea9aabab4ffc013811df72a76269d65fa Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 25 Feb 2020 20:39:19 +0000 Subject: [PATCH 05/10] Updates to plumes --- plumes.f90 | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/plumes.f90 b/plumes.f90 index 4714870..80b301f 100644 --- a/plumes.f90 +++ b/plumes.f90 @@ -1,24 +1,17 @@ subroutine plumes(V,L,AG,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. +!!January 2018 + +! The routine identifies all four-connected elements in a matrix. +! The returns are: +! V - arbitrary island number +! L - nx3 matrix where earch row has information about a particular cluster. +! the first column is the cluster number which corresponds to V. The second +! column is the number of elements in the cluster, and the third column is +! the value of the elements in that cluster (the elements of the input field). +! AG - Binary matrix with ones in the positions of the elements of the largers +! four-connected clusters. integer, intent(in) :: row,col,kend integer, intent(in), dimension(row,col) :: a From a52979cfffaf34b33c3244828bf4ceec281d5bd8 Mon Sep 17 00:00:00 2001 From: Philip Pegion Date: Tue, 17 Mar 2020 21:49:20 +0000 Subject: [PATCH 06/10] fixes for global_lats_h erroros --- atmosphere_stub.F90 | 289 ----------------------------------------- compile_standalone | 20 ++- setlats_lag_stochy.f | 29 +++-- standalone_stochy.F90 | 49 +++---- stochastic_physics.F90 | 1 + 5 files changed, 59 insertions(+), 329 deletions(-) diff --git a/atmosphere_stub.F90 b/atmosphere_stub.F90 index 38a7fd7..142504b 100644 --- a/atmosphere_stub.F90 +++ b/atmosphere_stub.F90 @@ -44,7 +44,6 @@ module atmosphere_stub_mod public :: atmosphere_init_stub !--- utility routines -!public :: atmosphere_return_winds, atmosphere_smooth_noise public :: atmosphere_resolution,atmosphere_domain,& atmosphere_scalar_field_halo,atmosphere_control_data @@ -175,111 +174,6 @@ subroutine atmosphere_init_stub (Grid_box, area) end subroutine atmosphere_init_stub -! subroutine atmosphere_smooth_noise (wnoise,npass,ns_type,renorm_type) -! -! !--- interface variables --- -! real,intent(inout) :: wnoise(isd:ied,jsd:jed,1) -! integer, intent(in) :: npass,ns_type,renorm_type -! !--- local variables -! integer:: i,j,nloops,nlast -! real ::inflation(isc:iec,jsc:jec),inflation2 -! ! scale factor for restoring inflation -! ! logic: -! ! if box mean: scalar get basic scaling, vector gets 1/grid dependent scaling 0-0 ; 0 - 1 -! ! if box mean2: no scaling -! ! if del2 : scalar gets grid dependent scaling,vector get basic scaling 1 0; 1 1 -! if(npass.GT.0) then -! if (ns_type.NE.2) then -! if (ns_type.EQ. 0) then -! !inflation2=1.0/sqrt(1.0/(4.0*npass)) -! inflation2=1.0/sqrt(1.0/(9.0*npass)) -! else -! inflation2=1.0/sqrt(1.0/(11.0/3.0*npass)) -! endif -! if ( ns_type.EQ.1) then ! del2 smoothing needs to be scaled by grid-size -! do j=jsc,jec -! do i=isc,iec -! inflation(i,j)=inflation2*Atm(mytile)%gridstruct%dxAV/(0.5*(Atm(mytile)%gridstruct%dx(i,j)+Atm(mytile)%gridstruct%dy(i,j))) -! enddo -! enddo -! else -! if ( renorm_type.EQ.1) then ! box smooth does not need scaling for scalar -! do j=jsc,jec -! do i=isc,iec -! inflation(i,j)=inflation2 -! enddo -! enddo -! else -! ! box mean needs inversize grid-size scaling for vector -! do j=jsc,jec -! do i=isc,iec -! inflation(i,j)=inflation2*(0.5*(Atm(mytile)%gridstruct%dx(i,j)+Atm(mytile)%gridstruct%dy(i,j)))/Atm(mytile)%gridstruct%dxAV -! enddo -! enddo -! endif -! endif -! endif -! nloops=npass/3 -! nlast=mod(npass,3) -! do j=1,nloops -! if (ns_type.EQ.1) then -! !call del2_cubed(wnoise , 0.25*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & -! call del2_cubed(wnoise , 0.20*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & -! Atm(mytile)%domain, npx, npy, 1, 3, Atm(mytile)%bd) -! else if (ns_type .EQ. 0) then -! call box_mean(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, 3, Atm(mytile)%bd) -! else if (ns_type .EQ. 2) then -! call box_mean2(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, 3, Atm(mytile)%bd) -! endif -! enddo -! if(nlast>0) then -! if (ns_type.EQ.1) then -! !call del2_cubed(wnoise , 0.25*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & -! call del2_cubed(wnoise , 0.20*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & -! Atm(mytile)%domain, npx, npy, 1, nlast, Atm(mytile)%bd) -! else if (ns_type .EQ. 0) then -! call box_mean(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, nlast, Atm(mytile)%bd) -! else if (ns_type .EQ. 2) then -! call box_mean2(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, nlast, Atm(mytile)%bd) -! endif -! endif -! ! restore amplitude -! if (ns_type.NE.2) then -! do j=jsc,jec -! do i=isc,iec -! wnoise(i,j,1)=wnoise(i,j,1)*inflation(i,j) -! enddo -! enddo -! endif -! endif -! end subroutine atmosphere_smooth_noise - -! subroutine atmosphere_return_winds (psi,ua,va,edge,km,vwts) -! integer,intent(in) :: edge -! real,intent(inout) :: psi(isd:ied,jsd:jed) -! real,intent(inout) :: ua(isc:iec+edge,jsc:jec) -! real,intent(inout) :: va(isc:iec,jsc:jec+edge) -! integer, optional,intent(in):: km -! real, optional,intent(in):: vwts(:) -! integer :: k -! call timing_on('COMM_TOTAL') -! call mpp_update_domains(psi, Atm(mytile)%domain, complete=.true.) -! call timing_off('COMM_TOTAL') -! if (edge.EQ.0) then -! call make_a_winds(ua, va, psi,Atm(mytile)%ng,Atm(mytile)%gridstruct,Atm(mytile)%bd,Atm(mytile)%npx,Atm(mytile)%npy) -! endif -! if (edge.EQ.1) then -! call make_c_winds(ua, va, psi,Atm(mytile)%ng,Atm(mytile)%gridstruct,Atm(mytile)%bd,Atm(mytile)%npx,Atm(mytile)%npy) -!! populate wind perturbations right here -! do k=1,km -! Atm(mytile)%urandom_c(isc:iec+edge,jsc:jec ,k)=ua*vwts(k) -! Atm(mytile)%vrandom_c(isc:iec ,jsc:jec+edge,k)=va*vwts(k) -! enddo -! !call mpp_update_domains(Atm(mytile)%urandom_c, Atm(mytile)%domain, complete=.true.) -! !call mpp_update_domains(Atm(mytile)%vrandom_c, Atm(mytile)%domain, complete=.true.) -! endif -! end subroutine atmosphere_return_winds -! subroutine del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd) !--------------------------------------------------------------- ! This routine is for filtering the omega field for the physics @@ -382,189 +276,6 @@ subroutine del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd) end subroutine del2_cubed -!>@brief The subroutine 'box_mean' filters a field with a 9-point mean stencil - - subroutine box_mean(q, gridstruct, domain, npx, npy, km, nmax, bd) - !--------------------------------------------------------------- - ! This routine is for filtering the omega field for the physics - !--------------------------------------------------------------- - integer, intent(in):: npx, npy, km, nmax - type(fv_grid_bounds_type), intent(IN) :: bd - real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km) - type(fv_grid_type), intent(IN), target :: gridstruct - type(domain2d), intent(INOUT) :: domain - real, parameter:: r3 = 1./3.,r9=1./9. - real :: q2(bd%isd:bd%ied,bd%jsd:bd%jed) - integer i,j,k, n, nt, ntimes - integer :: is, ie, js, je - integer :: isd, ied, jsd, jed - - !Local routine pointers -! real, pointer, dimension(:,:) :: rarea -! real, pointer, dimension(:,:) :: del6_u, del6_v -! logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner - - is = bd%is - ie = bd%ie - js = bd%js - je = bd%je - isd = bd%isd - ied = bd%ied - jsd = bd%jsd - jed = bd%jed - - ntimes = min(3, nmax) - - call timing_on('COMM_TOTAL') - call mpp_update_domains(q, domain, complete=.true.) - call timing_off('COMM_TOTAL') - - - do n=1,ntimes - nt = ntimes !- n - -!$OMP parallel do default(none) shared(km,is,ie,js,je,npx,npy, & -!$OMP q,nt,isd,jsd,gridstruct,bd) & -!$OMP private(q2) - do k=1,km - - if ( gridstruct%sw_corner ) then - q(1,1,k) = (q(1,1,k)+q(0,1,k)+q(1,0,k)) * r3 - q(0,1,k) = q(1,1,k) - q(1,0,k) = q(1,1,k) - endif - if ( gridstruct%se_corner ) then - q(ie, 1,k) = (q(ie,1,k)+q(npx,1,k)+q(ie,0,k)) * r3 - q(npx,1,k) = q(ie,1,k) - q(ie, 0,k) = q(ie,1,k) - endif - if ( gridstruct%ne_corner ) then - q(ie, je,k) = (q(ie,je,k)+q(npx,je,k)+q(ie,npy,k)) * r3 - q(npx,je,k) = q(ie,je,k) - q(ie,npy,k) = q(ie,je,k) - endif - if ( gridstruct%nw_corner ) then - q(1, je,k) = (q(1,je,k)+q(0,je,k)+q(1,npy,k)) * r3 - q(0, je,k) = q(1,je,k) - q(1,npy,k) = q(1,je,k) - endif - - if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1, gridstruct%nested, bd, & - gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner ) - - if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2, gridstruct%nested, bd, & - gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) - - !do j=js-nt,je+nt - ! do i=is-nt,ie+nt - do j=jsd+1,jed-1 - do i=isd+1,ied-1 - !q2(i,j) = (gridstruct%area(i-1,j+1)*q(i-1,j+1,k) + gridstruct%area(i,j+1)*q(i,j+1,k) + gridstruct%area(i+1,j+1)*q(i+1,j+1,k) +& - ! gridstruct%area(i-1,j )*q(i-1,j,k) + gridstruct%area(i,j )*q(i,j ,k) + gridstruct%area(i+1,j )*q(i+1,j ,k) +& - ! gridstruct%area(i-1,j-1)*q(i-1,j-1,k) + gridstruct%area(i,j-1)*q(i,j-1,k) + gridstruct%area(i+1,j-1)*q(i+1,j-1,k))/SUM(gridstruct%area(i-1:i+1,j-1:j+1)) - q2(i,j) = r9*(q(i-1,j+1,k)+q(i,j+1,k)+q(i+1,j+1,k)+q(i-1,j,k)+q(i,j,k)+q(i+1,j,k)+q(i-1,j-1,k)+q(i,j-1,k)+q(i+1,j-1,k)) - !if (j.GE. je .AND. i.GE. ie) print*,'area +1=',gridstruct%area(i-1:i+1,j+1) - !if (j.GE. je .AND. i.GE. ie) print*,'area =',gridstruct%area(i-1:i+1,j) - !if (j.GE. je .AND. i.GE. ie) print*,'area -1=',gridstruct%area(i-1:i+1,j-1) - !if (j.GE. je .AND. i.GE. ie) print*,'q +1=',q(i-1:i+1,j+1,k) - !if (j.GE. je .AND. i.GE. ie) print*,'q =',q(i-1:i+1,j,k) - !if (j.GE. je .AND. i.GE. ie) print*,'q -1=',q(i-1:i+1,j-1,k) - enddo - enddo - do j=js-nt,je+nt - do i=is-nt,ie+nt - q(i,j,k)=q2(i,j) - enddo - enddo - enddo - enddo - end subroutine box_mean - - subroutine box_mean2(q, gridstruct, domain, npx, npy, km, nmax, bd) - !--------------------------------------------------------------- - ! This routine is for filtering the omega field for the physics - !--------------------------------------------------------------- - integer, intent(in):: npx, npy, km, nmax - type(fv_grid_bounds_type), intent(IN) :: bd - real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km) - type(fv_grid_type), intent(IN), target :: gridstruct - type(domain2d), intent(INOUT) :: domain - real, parameter:: r3 = 1./3.,r10=0.1 - real :: q2(bd%isd:bd%ied,bd%jsd:bd%jed) - integer i,j,k, n, nt, ntimes - integer :: is, ie, js, je - integer :: isd, ied, jsd, jed - - !Local routine pointers -! real, pointer, dimension(:,:) :: rarea -! real, pointer, dimension(:,:) :: del6_u, del6_v -! logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner - - is = bd%is - ie = bd%ie - js = bd%js - je = bd%je - isd = bd%isd - ied = bd%ied - jsd = bd%jsd - jed = bd%jed - - ntimes = min(3, nmax) - - call timing_on('COMM_TOTAL') - call mpp_update_domains(q, domain, complete=.true.) - call timing_off('COMM_TOTAL') - - - do n=1,ntimes - nt = ntimes !- n - -!$OMP parallel do default(none) shared(km,is,ie,js,je,npx,npy, & -!$OMP q,nt,isd,jsd,gridstruct,bd) & -!$OMP private(q2) - do k=1,km - - if ( gridstruct%sw_corner ) then - q(1,1,k) = (q(1,1,k)+q(0,1,k)+q(1,0,k)) * r3 - q(0,1,k) = q(1,1,k) - q(1,0,k) = q(1,1,k) - endif - if ( gridstruct%se_corner ) then - q(ie, 1,k) = (q(ie,1,k)+q(npx,1,k)+q(ie,0,k)) * r3 - q(npx,1,k) = q(ie,1,k) - q(ie, 0,k) = q(ie,1,k) - endif - if ( gridstruct%ne_corner ) then - q(ie, je,k) = (q(ie,je,k)+q(npx,je,k)+q(ie,npy,k)) * r3 - q(npx,je,k) = q(ie,je,k) - q(ie,npy,k) = q(ie,je,k) - endif - if ( gridstruct%nw_corner ) then - q(1, je,k) = (q(1,je,k)+q(0,je,k)+q(1,npy,k)) * r3 - q(0, je,k) = q(1,je,k) - q(1,npy,k) = q(1,je,k) - endif - - if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1, gridstruct%nested, bd, & - gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner ) - - if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2, gridstruct%nested, bd, & - gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) - - do j=jsd+1,jed-1 - do i=isd+1,ied-1 - q2(i,j) = r10*(q(i-1,j+1,k)+q(i,j+1,k)+q(i+1,j+1,k)+q(i-1,j,k)+2*q(i,j,k)+q(i+1,j,k)+q(i-1,j-1,k)+q(i,j-1,k)+q(i+1,j-1,k)) - enddo - enddo - do j=js-nt,je+nt - do i=is-nt,ie+nt - q(i,j,k)=q2(i,j) - enddo - enddo - enddo - enddo - - end subroutine box_mean2 subroutine make_a_winds(ua, va, psi, ng, gridstruct, bd, npx, npy) integer, intent(IN) :: ng, npx, npy diff --git a/compile_standalone b/compile_standalone index 11daf3e..7598a6e 100755 --- a/compile_standalone +++ b/compile_standalone @@ -1,13 +1,15 @@ #!/bin/sh -x -compile_all=1 +compile_all=0 +SRC_DIR=/scratch2/BMC/gsienkf/Philip.Pegion/ufs-s2s-model/ FC=mpif90 -INCS="-I. -I../FV3/gfsphysics/ -I../FMS/include -I../FV3/atmos_cubed_sphere -I../FMS/fv3gfs" -FLAGS64=" -traceback -real-size 64 -DSTOCHY_UNIT_TEST -c "$INCS -FLAGS=" -traceback -DSTOCHY_UNIT_TEST -c "$INCS +OPT_FLAGS="-g -O0 -check all -check noarg_temp_created -check nopointer -warn -warn noerrors -fp-stack-check -fstack-protector-all -fpe0 -debug -traceback -ftrapuv" +INCS="-I. -I${SRC_DIR}FV3/gfsphysics/ -I${SRC_DIR}/FMS/include -I${SRC_DIR}/FV3/atmos_cubed_sphere -I${SRC_DIR}/FMS/fv3gfs" +FLAGS64=$OPT_FLAGS" -qopenmp -real-size 64 -DSTOCHY_UNIT_TEST -c "$INCS +FLAGS=$OPT_FLAGS" -qopenmp -DSTOCHY_UNIT_TEST -c "$INCS if [ $compile_all -eq 1 ];then rm -f *.i90 *.i *.o *.mod lib*a - $FC ${FLAGS} fv_control_stub.F90 - $FC ${FLAGS} atmosphere_stub.F90 + $FC ${FLAGS64} fv_control_stub.F90 + $FC ${FLAGS64} atmosphere_stub.F90 $FC ${FLAGS64} stochy_namelist_def.F90 $FC ${FLAGS64} stochy_resol_def.f $FC ${FLAGS64} stochy_gg_def.f @@ -38,4 +40,8 @@ if [ $compile_all -eq 1 ];then $FC ${FLAGS64} stochastic_physics.F90 ar rv libstochastic_physics.a *.o fi -$FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf +rm libstochastic_physics.a +$FC ${FLAGS64} setlats_lag_stochy.f +$FC ${FLAGS64} stochastic_physics.F90 +ar rv libstochastic_physics.a *.o +$FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L${SRC_DIR}/FV3/atmos_cubed_sphere -lfv3core -L${SRC_DIR}/FMS/FMS_INSTALL -lfms -L${SRC_DIR}/FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf diff --git a/setlats_lag_stochy.f b/setlats_lag_stochy.f index ff4e4f0..98f70e3 100644 --- a/setlats_lag_stochy.f +++ b/setlats_lag_stochy.f @@ -54,26 +54,35 @@ subroutine setlats_lag_stochy(lats_nodes_a, global_lats_a, global_lats_h(j1-j2) = global_lats_a(latg) + 1 - jj ! set south pole yhalo enddo ! + !if (me==0) print*,'lats_nodes_a=',lats_nodes_a(1),latg + !if (me==0) print*,'lats_nodes_h=',lats_nodes_h(:) + if (lats_nodes_a(1) /= latg) then ! ! set non-polar south yhalos jpt_h = 0 do nn=1,nodes-1 - jpt_h = jpt_h + lats_nodes_h(nn) - lat_val = global_lats_h(jpt_h-yhalo) - do jj=1,yhalo - global_lats_h(jpt_h-yhalo+jj) = min(lat_val+jj,latg) - enddo + if(lats_nodes_h(nn).GT.0) then + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h-yhalo) + do jj=1,yhalo + global_lats_h(jpt_h-yhalo+jj) = min(lat_val+jj,latg) + enddo + endif enddo ! ! set non-polar north yhalos jpt_h = 0 + !if (me==0) print*,'in setlats',nodes,yhalo,size(lats_nodes_h) do nn=1,nodes-1 - jpt_h = jpt_h + lats_nodes_h(nn) - lat_val = global_lats_h(jpt_h+yhalo+1) - do jj=1,yhalo - global_lats_h(jpt_h+yhalo-(jj-1)) = max(lat_val-jj,1) - enddo + ! if (me==0) print*,'in loop',nn,jpt_h,lats_nodes_h(nn) + if(lats_nodes_h(nn).GT.0) then + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h+yhalo+1) + do jj=1,yhalo + global_lats_h(jpt_h+yhalo-(jj-1)) = max(lat_val-jj,1) + enddo + endif enddo ! endif diff --git a/standalone_stochy.F90 b/standalone_stochy.F90 index 75e9ce7..67de8e2 100644 --- a/standalone_stochy.F90 +++ b/standalone_stochy.F90 @@ -23,7 +23,7 @@ program standalone_stochy integer :: ncid,xt_dim_id,yt_dim_id,time_dim_id,xt_var_id,yt_var_id,time_var_id integer :: varid1,varid2,varid3,varid4 integer :: zt_dim_id,zt_var_id -character*1 :: strid +character*3 :: strid type(GFS_grid_type),allocatable :: Grid(:) type(GFS_coupling_type),allocatable :: Coupling(:) ! stochastic namelist fields @@ -37,9 +37,10 @@ program standalone_stochy real(kind=kind_dbl_prec) fhstoch,skeb_diss_smooth,spptint,shumint,skebint,skebnorm real(kind=kind_dbl_prec), dimension(5) :: skeb,skeb_lscale,skeb_tau real(kind=kind_dbl_prec), dimension(5) :: sppt,sppt_lscale,sppt_tau +real(kind=kind_dbl_prec), dimension(5) :: ocnp,ocnp_lscale,ocnp_tau real(kind=kind_dbl_prec), dimension(5) :: shum,shum_lscale,shum_tau integer,dimension(5) ::skeb_vfilt -integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb +integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_ocnp logical stochini,sppt_logit,new_lscale logical use_zmtnblck logical sppt_land @@ -73,7 +74,7 @@ program standalone_stochy logical :: write_this_tile integer :: nargs,ntile_out,nlunit,pe,npes,stackmax=4000000 character*80 :: fname -character*1 :: ntile_out_str +character*2 :: ntile_out_str real(kind=4),allocatable,dimension(:,:) :: workg real(kind=4),allocatable,dimension(:,:,:) :: workg3d @@ -88,14 +89,15 @@ program standalone_stochy 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,& - shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale + shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & + ocnp,ocnp_lscale,ocnp_tau,iseed_ocnp write_this_tile=.false. ntile_out_str='0' nargs=iargc() if (nargs.EQ.1) then call getarg(1,ntile_out_str) endif -read(ntile_out_str,'(I1.1)') ntile_out +read(ntile_out_str,'(I2.2)') ntile_out open (unit=nlunit, file='input.nml', READONLY, status='OLD') read(nlunit,nam_stochy) close(nlunit) @@ -139,10 +141,12 @@ program standalone_stochy jec=Atm(1)%bd%jec nx=Atm(1)%npx-1 ny=Atm(1)%npy-1 -allocate(workg(nx,ny)) -allocate(workg3d(nx,ny,nlevs)) -nblks=ny -blksz=nx +nx2=iec-isc+1 +ny2=jec-jsc+1 +allocate(workg(nx2,ny2)) +allocate(workg3d(nx2,ny2,nlevs)) +blksz=nx2 +nblks=ny2 Allocate(Model%blksz(nblks)) Model%blksz(:)=blksz nthreads = omp_get_num_threads() @@ -180,27 +184,26 @@ program standalone_stochy Grid(j)%xlat(:)=Init_parm%xlat(:,j) Grid(j)%xlon(:)=Init_parm%xlon(:,j) enddo -allocate(grid_xt(nx),grid_yt(ny)) -do i=1,nx +allocate(grid_xt(nx2),grid_yt(ny2)) +do i=1,nx2 grid_xt(i)=i enddo -do i=1,ny +do i=1,ny2 grid_yt(i)=i enddo !setup GFS_coupling allocate(Coupling(nblks)) call init_stochastic_physics(Model, Init_parm, ntasks, nthreads) call get_outfile(fname) -write(strid,'(I1.1)') my_id+1 +write(strid,'(I2.1)') my_id+1 if (ntile_out.EQ.0) write_this_tile=.true. if ((my_id+1).EQ.ntile_out) write_this_tile=.true. -print*,trim(fname)//'.tile'//strid//'.nc',write_this_tile +print*,trim(fname)//'.tile'//trim(adjustl(strid))//'.nc',write_this_tile if (write_this_tile) then fid=30+my_id -!ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) -ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) -ierr=NF90_DEF_DIM(ncid,"grid_xt",nx,xt_dim_id) -ierr=NF90_DEF_DIM(ncid,"grid_yt",ny,yt_dim_id) +ierr=nf90_create(trim(fname)//'.tile'//trim(adjustl(strid))//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +ierr=NF90_DEF_DIM(ncid,"grid_xt",nx2,xt_dim_id) +ierr=NF90_DEF_DIM(ncid,"grid_yt",ny2,yt_dim_id) if (Model%do_skeb)ierr=NF90_DEF_DIM(ncid,"p_ref",nlevs,zt_dim_id) ierr=NF90_DEF_DIM(ncid,"time",NF90_UNLIMITED,time_dim_id) !> - Define the dimension variables. @@ -269,33 +272,33 @@ program standalone_stochy if (Model%do_skeb)allocate(Coupling(i)%skebu_wts(blksz,nlevs)) if (Model%do_skeb)allocate(Coupling(i)%skebv_wts(blksz,nlevs)) enddo -do i=1,200 +do i=1,10 Model%kdt=i ts=i/4.0 call run_stochastic_physics(Model, Grid, Coupling, nthreads) if (Model%me.EQ.0) print*,'sppt_wts=',i,Coupling(1)%sppt_wts(1,20) if (write_this_tile) then if (Model%do_sppt)then - do j=1,ny + do j=1,ny2 workg(:,j)=Coupling(j)%sppt_wts(:,20) enddo ierr=NF90_PUT_VAR(ncid,varid1,workg,(/1,1,i/)) endif if (Model%do_shum)then - do j=1,ny + do j=1,ny2 workg(:,j)=Coupling(j)%shum_wts(:,1) enddo ierr=NF90_PUT_VAR(ncid,varid2,workg,(/1,1,i/)) endif if (Model%do_skeb)then do k=1,nlevs - do j=1,ny + do j=1,ny2 workg3d(:,j,k)=Coupling(j)%skebu_wts(:,k) enddo enddo ierr=NF90_PUT_VAR(ncid,varid3,workg3d,(/1,1,1,i/)) do k=1,nlevs - do j=1,ny + do j=1,ny2 workg3d(:,j,k)=Coupling(j)%skebv_wts(:,k) enddo enddo diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 9b89902..a09e5d9 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -59,6 +59,7 @@ subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) nodes=ntasks gis_stochy%me=me gis_stochy%nodes=nodes +gis_stochy%yhalo=10 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 From 21cc53974a94d6ba7bb336fb72b820f80a15a12a Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 25 Mar 2020 15:30:24 +0000 Subject: [PATCH 07/10] bringing back these routines to psd/develop --- atmosphere_stub.F90 | 289 +++++++++++++++++++++++++++++++++++++++++ compile_standalone | 20 +-- setlats_lag_stochy.f | 29 ++--- standalone_stochy.F90 | 49 ++++--- stochastic_physics.F90 | 1 - 5 files changed, 329 insertions(+), 59 deletions(-) diff --git a/atmosphere_stub.F90 b/atmosphere_stub.F90 index 142504b..38a7fd7 100644 --- a/atmosphere_stub.F90 +++ b/atmosphere_stub.F90 @@ -44,6 +44,7 @@ module atmosphere_stub_mod public :: atmosphere_init_stub !--- utility routines +!public :: atmosphere_return_winds, atmosphere_smooth_noise public :: atmosphere_resolution,atmosphere_domain,& atmosphere_scalar_field_halo,atmosphere_control_data @@ -174,6 +175,111 @@ subroutine atmosphere_init_stub (Grid_box, area) end subroutine atmosphere_init_stub +! subroutine atmosphere_smooth_noise (wnoise,npass,ns_type,renorm_type) +! +! !--- interface variables --- +! real,intent(inout) :: wnoise(isd:ied,jsd:jed,1) +! integer, intent(in) :: npass,ns_type,renorm_type +! !--- local variables +! integer:: i,j,nloops,nlast +! real ::inflation(isc:iec,jsc:jec),inflation2 +! ! scale factor for restoring inflation +! ! logic: +! ! if box mean: scalar get basic scaling, vector gets 1/grid dependent scaling 0-0 ; 0 - 1 +! ! if box mean2: no scaling +! ! if del2 : scalar gets grid dependent scaling,vector get basic scaling 1 0; 1 1 +! if(npass.GT.0) then +! if (ns_type.NE.2) then +! if (ns_type.EQ. 0) then +! !inflation2=1.0/sqrt(1.0/(4.0*npass)) +! inflation2=1.0/sqrt(1.0/(9.0*npass)) +! else +! inflation2=1.0/sqrt(1.0/(11.0/3.0*npass)) +! endif +! if ( ns_type.EQ.1) then ! del2 smoothing needs to be scaled by grid-size +! do j=jsc,jec +! do i=isc,iec +! inflation(i,j)=inflation2*Atm(mytile)%gridstruct%dxAV/(0.5*(Atm(mytile)%gridstruct%dx(i,j)+Atm(mytile)%gridstruct%dy(i,j))) +! enddo +! enddo +! else +! if ( renorm_type.EQ.1) then ! box smooth does not need scaling for scalar +! do j=jsc,jec +! do i=isc,iec +! inflation(i,j)=inflation2 +! enddo +! enddo +! else +! ! box mean needs inversize grid-size scaling for vector +! do j=jsc,jec +! do i=isc,iec +! inflation(i,j)=inflation2*(0.5*(Atm(mytile)%gridstruct%dx(i,j)+Atm(mytile)%gridstruct%dy(i,j)))/Atm(mytile)%gridstruct%dxAV +! enddo +! enddo +! endif +! endif +! endif +! nloops=npass/3 +! nlast=mod(npass,3) +! do j=1,nloops +! if (ns_type.EQ.1) then +! !call del2_cubed(wnoise , 0.25*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & +! call del2_cubed(wnoise , 0.20*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & +! Atm(mytile)%domain, npx, npy, 1, 3, Atm(mytile)%bd) +! else if (ns_type .EQ. 0) then +! call box_mean(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, 3, Atm(mytile)%bd) +! else if (ns_type .EQ. 2) then +! call box_mean2(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, 3, Atm(mytile)%bd) +! endif +! enddo +! if(nlast>0) then +! if (ns_type.EQ.1) then +! !call del2_cubed(wnoise , 0.25*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & +! call del2_cubed(wnoise , 0.20*Atm(mytile)%gridstruct%da_min, Atm(mytile)%gridstruct, & +! Atm(mytile)%domain, npx, npy, 1, nlast, Atm(mytile)%bd) +! else if (ns_type .EQ. 0) then +! call box_mean(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, nlast, Atm(mytile)%bd) +! else if (ns_type .EQ. 2) then +! call box_mean2(wnoise , Atm(mytile)%gridstruct, Atm(mytile)%domain, Atm(mytile)%npx, Atm(mytile)%npy, 1, nlast, Atm(mytile)%bd) +! endif +! endif +! ! restore amplitude +! if (ns_type.NE.2) then +! do j=jsc,jec +! do i=isc,iec +! wnoise(i,j,1)=wnoise(i,j,1)*inflation(i,j) +! enddo +! enddo +! endif +! endif +! end subroutine atmosphere_smooth_noise + +! subroutine atmosphere_return_winds (psi,ua,va,edge,km,vwts) +! integer,intent(in) :: edge +! real,intent(inout) :: psi(isd:ied,jsd:jed) +! real,intent(inout) :: ua(isc:iec+edge,jsc:jec) +! real,intent(inout) :: va(isc:iec,jsc:jec+edge) +! integer, optional,intent(in):: km +! real, optional,intent(in):: vwts(:) +! integer :: k +! call timing_on('COMM_TOTAL') +! call mpp_update_domains(psi, Atm(mytile)%domain, complete=.true.) +! call timing_off('COMM_TOTAL') +! if (edge.EQ.0) then +! call make_a_winds(ua, va, psi,Atm(mytile)%ng,Atm(mytile)%gridstruct,Atm(mytile)%bd,Atm(mytile)%npx,Atm(mytile)%npy) +! endif +! if (edge.EQ.1) then +! call make_c_winds(ua, va, psi,Atm(mytile)%ng,Atm(mytile)%gridstruct,Atm(mytile)%bd,Atm(mytile)%npx,Atm(mytile)%npy) +!! populate wind perturbations right here +! do k=1,km +! Atm(mytile)%urandom_c(isc:iec+edge,jsc:jec ,k)=ua*vwts(k) +! Atm(mytile)%vrandom_c(isc:iec ,jsc:jec+edge,k)=va*vwts(k) +! enddo +! !call mpp_update_domains(Atm(mytile)%urandom_c, Atm(mytile)%domain, complete=.true.) +! !call mpp_update_domains(Atm(mytile)%vrandom_c, Atm(mytile)%domain, complete=.true.) +! endif +! end subroutine atmosphere_return_winds +! subroutine del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd) !--------------------------------------------------------------- ! This routine is for filtering the omega field for the physics @@ -276,6 +382,189 @@ subroutine del2_cubed(q, cd, gridstruct, domain, npx, npy, km, nmax, bd) end subroutine del2_cubed +!>@brief The subroutine 'box_mean' filters a field with a 9-point mean stencil + + subroutine box_mean(q, gridstruct, domain, npx, npy, km, nmax, bd) + !--------------------------------------------------------------- + ! This routine is for filtering the omega field for the physics + !--------------------------------------------------------------- + integer, intent(in):: npx, npy, km, nmax + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km) + type(fv_grid_type), intent(IN), target :: gridstruct + type(domain2d), intent(INOUT) :: domain + real, parameter:: r3 = 1./3.,r9=1./9. + real :: q2(bd%isd:bd%ied,bd%jsd:bd%jed) + integer i,j,k, n, nt, ntimes + integer :: is, ie, js, je + integer :: isd, ied, jsd, jed + + !Local routine pointers +! real, pointer, dimension(:,:) :: rarea +! real, pointer, dimension(:,:) :: del6_u, del6_v +! logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + + ntimes = min(3, nmax) + + call timing_on('COMM_TOTAL') + call mpp_update_domains(q, domain, complete=.true.) + call timing_off('COMM_TOTAL') + + + do n=1,ntimes + nt = ntimes !- n + +!$OMP parallel do default(none) shared(km,is,ie,js,je,npx,npy, & +!$OMP q,nt,isd,jsd,gridstruct,bd) & +!$OMP private(q2) + do k=1,km + + if ( gridstruct%sw_corner ) then + q(1,1,k) = (q(1,1,k)+q(0,1,k)+q(1,0,k)) * r3 + q(0,1,k) = q(1,1,k) + q(1,0,k) = q(1,1,k) + endif + if ( gridstruct%se_corner ) then + q(ie, 1,k) = (q(ie,1,k)+q(npx,1,k)+q(ie,0,k)) * r3 + q(npx,1,k) = q(ie,1,k) + q(ie, 0,k) = q(ie,1,k) + endif + if ( gridstruct%ne_corner ) then + q(ie, je,k) = (q(ie,je,k)+q(npx,je,k)+q(ie,npy,k)) * r3 + q(npx,je,k) = q(ie,je,k) + q(ie,npy,k) = q(ie,je,k) + endif + if ( gridstruct%nw_corner ) then + q(1, je,k) = (q(1,je,k)+q(0,je,k)+q(1,npy,k)) * r3 + q(0, je,k) = q(1,je,k) + q(1,npy,k) = q(1,je,k) + endif + + if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner ) + + if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) + + !do j=js-nt,je+nt + ! do i=is-nt,ie+nt + do j=jsd+1,jed-1 + do i=isd+1,ied-1 + !q2(i,j) = (gridstruct%area(i-1,j+1)*q(i-1,j+1,k) + gridstruct%area(i,j+1)*q(i,j+1,k) + gridstruct%area(i+1,j+1)*q(i+1,j+1,k) +& + ! gridstruct%area(i-1,j )*q(i-1,j,k) + gridstruct%area(i,j )*q(i,j ,k) + gridstruct%area(i+1,j )*q(i+1,j ,k) +& + ! gridstruct%area(i-1,j-1)*q(i-1,j-1,k) + gridstruct%area(i,j-1)*q(i,j-1,k) + gridstruct%area(i+1,j-1)*q(i+1,j-1,k))/SUM(gridstruct%area(i-1:i+1,j-1:j+1)) + q2(i,j) = r9*(q(i-1,j+1,k)+q(i,j+1,k)+q(i+1,j+1,k)+q(i-1,j,k)+q(i,j,k)+q(i+1,j,k)+q(i-1,j-1,k)+q(i,j-1,k)+q(i+1,j-1,k)) + !if (j.GE. je .AND. i.GE. ie) print*,'area +1=',gridstruct%area(i-1:i+1,j+1) + !if (j.GE. je .AND. i.GE. ie) print*,'area =',gridstruct%area(i-1:i+1,j) + !if (j.GE. je .AND. i.GE. ie) print*,'area -1=',gridstruct%area(i-1:i+1,j-1) + !if (j.GE. je .AND. i.GE. ie) print*,'q +1=',q(i-1:i+1,j+1,k) + !if (j.GE. je .AND. i.GE. ie) print*,'q =',q(i-1:i+1,j,k) + !if (j.GE. je .AND. i.GE. ie) print*,'q -1=',q(i-1:i+1,j-1,k) + enddo + enddo + do j=js-nt,je+nt + do i=is-nt,ie+nt + q(i,j,k)=q2(i,j) + enddo + enddo + enddo + enddo + end subroutine box_mean + + subroutine box_mean2(q, gridstruct, domain, npx, npy, km, nmax, bd) + !--------------------------------------------------------------- + ! This routine is for filtering the omega field for the physics + !--------------------------------------------------------------- + integer, intent(in):: npx, npy, km, nmax + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed,km) + type(fv_grid_type), intent(IN), target :: gridstruct + type(domain2d), intent(INOUT) :: domain + real, parameter:: r3 = 1./3.,r10=0.1 + real :: q2(bd%isd:bd%ied,bd%jsd:bd%jed) + integer i,j,k, n, nt, ntimes + integer :: is, ie, js, je + integer :: isd, ied, jsd, jed + + !Local routine pointers +! real, pointer, dimension(:,:) :: rarea +! real, pointer, dimension(:,:) :: del6_u, del6_v +! logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + + ntimes = min(3, nmax) + + call timing_on('COMM_TOTAL') + call mpp_update_domains(q, domain, complete=.true.) + call timing_off('COMM_TOTAL') + + + do n=1,ntimes + nt = ntimes !- n + +!$OMP parallel do default(none) shared(km,is,ie,js,je,npx,npy, & +!$OMP q,nt,isd,jsd,gridstruct,bd) & +!$OMP private(q2) + do k=1,km + + if ( gridstruct%sw_corner ) then + q(1,1,k) = (q(1,1,k)+q(0,1,k)+q(1,0,k)) * r3 + q(0,1,k) = q(1,1,k) + q(1,0,k) = q(1,1,k) + endif + if ( gridstruct%se_corner ) then + q(ie, 1,k) = (q(ie,1,k)+q(npx,1,k)+q(ie,0,k)) * r3 + q(npx,1,k) = q(ie,1,k) + q(ie, 0,k) = q(ie,1,k) + endif + if ( gridstruct%ne_corner ) then + q(ie, je,k) = (q(ie,je,k)+q(npx,je,k)+q(ie,npy,k)) * r3 + q(npx,je,k) = q(ie,je,k) + q(ie,npy,k) = q(ie,je,k) + endif + if ( gridstruct%nw_corner ) then + q(1, je,k) = (q(1,je,k)+q(0,je,k)+q(1,npy,k)) * r3 + q(0, je,k) = q(1,je,k) + q(1,npy,k) = q(1,je,k) + endif + + if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner ) + + if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2, gridstruct%nested, bd, & + gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) + + do j=jsd+1,jed-1 + do i=isd+1,ied-1 + q2(i,j) = r10*(q(i-1,j+1,k)+q(i,j+1,k)+q(i+1,j+1,k)+q(i-1,j,k)+2*q(i,j,k)+q(i+1,j,k)+q(i-1,j-1,k)+q(i,j-1,k)+q(i+1,j-1,k)) + enddo + enddo + do j=js-nt,je+nt + do i=is-nt,ie+nt + q(i,j,k)=q2(i,j) + enddo + enddo + enddo + enddo + + end subroutine box_mean2 subroutine make_a_winds(ua, va, psi, ng, gridstruct, bd, npx, npy) integer, intent(IN) :: ng, npx, npy diff --git a/compile_standalone b/compile_standalone index 7598a6e..11daf3e 100755 --- a/compile_standalone +++ b/compile_standalone @@ -1,15 +1,13 @@ #!/bin/sh -x -compile_all=0 -SRC_DIR=/scratch2/BMC/gsienkf/Philip.Pegion/ufs-s2s-model/ +compile_all=1 FC=mpif90 -OPT_FLAGS="-g -O0 -check all -check noarg_temp_created -check nopointer -warn -warn noerrors -fp-stack-check -fstack-protector-all -fpe0 -debug -traceback -ftrapuv" -INCS="-I. -I${SRC_DIR}FV3/gfsphysics/ -I${SRC_DIR}/FMS/include -I${SRC_DIR}/FV3/atmos_cubed_sphere -I${SRC_DIR}/FMS/fv3gfs" -FLAGS64=$OPT_FLAGS" -qopenmp -real-size 64 -DSTOCHY_UNIT_TEST -c "$INCS -FLAGS=$OPT_FLAGS" -qopenmp -DSTOCHY_UNIT_TEST -c "$INCS +INCS="-I. -I../FV3/gfsphysics/ -I../FMS/include -I../FV3/atmos_cubed_sphere -I../FMS/fv3gfs" +FLAGS64=" -traceback -real-size 64 -DSTOCHY_UNIT_TEST -c "$INCS +FLAGS=" -traceback -DSTOCHY_UNIT_TEST -c "$INCS if [ $compile_all -eq 1 ];then rm -f *.i90 *.i *.o *.mod lib*a - $FC ${FLAGS64} fv_control_stub.F90 - $FC ${FLAGS64} atmosphere_stub.F90 + $FC ${FLAGS} fv_control_stub.F90 + $FC ${FLAGS} atmosphere_stub.F90 $FC ${FLAGS64} stochy_namelist_def.F90 $FC ${FLAGS64} stochy_resol_def.f $FC ${FLAGS64} stochy_gg_def.f @@ -40,8 +38,4 @@ if [ $compile_all -eq 1 ];then $FC ${FLAGS64} stochastic_physics.F90 ar rv libstochastic_physics.a *.o fi -rm libstochastic_physics.a -$FC ${FLAGS64} setlats_lag_stochy.f -$FC ${FLAGS64} stochastic_physics.F90 -ar rv libstochastic_physics.a *.o -$FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L${SRC_DIR}/FV3/atmos_cubed_sphere -lfv3core -L${SRC_DIR}/FMS/FMS_INSTALL -lfms -L${SRC_DIR}/FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf +$FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf diff --git a/setlats_lag_stochy.f b/setlats_lag_stochy.f index 98f70e3..ff4e4f0 100644 --- a/setlats_lag_stochy.f +++ b/setlats_lag_stochy.f @@ -54,35 +54,26 @@ subroutine setlats_lag_stochy(lats_nodes_a, global_lats_a, global_lats_h(j1-j2) = global_lats_a(latg) + 1 - jj ! set south pole yhalo enddo ! - !if (me==0) print*,'lats_nodes_a=',lats_nodes_a(1),latg - !if (me==0) print*,'lats_nodes_h=',lats_nodes_h(:) - if (lats_nodes_a(1) /= latg) then ! ! set non-polar south yhalos jpt_h = 0 do nn=1,nodes-1 - if(lats_nodes_h(nn).GT.0) then - jpt_h = jpt_h + lats_nodes_h(nn) - lat_val = global_lats_h(jpt_h-yhalo) - do jj=1,yhalo - global_lats_h(jpt_h-yhalo+jj) = min(lat_val+jj,latg) - enddo - endif + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h-yhalo) + do jj=1,yhalo + global_lats_h(jpt_h-yhalo+jj) = min(lat_val+jj,latg) + enddo enddo ! ! set non-polar north yhalos jpt_h = 0 - !if (me==0) print*,'in setlats',nodes,yhalo,size(lats_nodes_h) do nn=1,nodes-1 - ! if (me==0) print*,'in loop',nn,jpt_h,lats_nodes_h(nn) - if(lats_nodes_h(nn).GT.0) then - jpt_h = jpt_h + lats_nodes_h(nn) - lat_val = global_lats_h(jpt_h+yhalo+1) - do jj=1,yhalo - global_lats_h(jpt_h+yhalo-(jj-1)) = max(lat_val-jj,1) - enddo - endif + jpt_h = jpt_h + lats_nodes_h(nn) + lat_val = global_lats_h(jpt_h+yhalo+1) + do jj=1,yhalo + global_lats_h(jpt_h+yhalo-(jj-1)) = max(lat_val-jj,1) + enddo enddo ! endif diff --git a/standalone_stochy.F90 b/standalone_stochy.F90 index 67de8e2..75e9ce7 100644 --- a/standalone_stochy.F90 +++ b/standalone_stochy.F90 @@ -23,7 +23,7 @@ program standalone_stochy integer :: ncid,xt_dim_id,yt_dim_id,time_dim_id,xt_var_id,yt_var_id,time_var_id integer :: varid1,varid2,varid3,varid4 integer :: zt_dim_id,zt_var_id -character*3 :: strid +character*1 :: strid type(GFS_grid_type),allocatable :: Grid(:) type(GFS_coupling_type),allocatable :: Coupling(:) ! stochastic namelist fields @@ -37,10 +37,9 @@ program standalone_stochy real(kind=kind_dbl_prec) fhstoch,skeb_diss_smooth,spptint,shumint,skebint,skebnorm real(kind=kind_dbl_prec), dimension(5) :: skeb,skeb_lscale,skeb_tau real(kind=kind_dbl_prec), dimension(5) :: sppt,sppt_lscale,sppt_tau -real(kind=kind_dbl_prec), dimension(5) :: ocnp,ocnp_lscale,ocnp_tau real(kind=kind_dbl_prec), dimension(5) :: shum,shum_lscale,shum_tau integer,dimension(5) ::skeb_vfilt -integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_ocnp +integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb logical stochini,sppt_logit,new_lscale logical use_zmtnblck logical sppt_land @@ -74,7 +73,7 @@ program standalone_stochy logical :: write_this_tile integer :: nargs,ntile_out,nlunit,pe,npes,stackmax=4000000 character*80 :: fname -character*2 :: ntile_out_str +character*1 :: ntile_out_str real(kind=4),allocatable,dimension(:,:) :: workg real(kind=4),allocatable,dimension(:,:,:) :: workg3d @@ -89,15 +88,14 @@ program standalone_stochy 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,& - shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & - ocnp,ocnp_lscale,ocnp_tau,iseed_ocnp + shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale write_this_tile=.false. ntile_out_str='0' nargs=iargc() if (nargs.EQ.1) then call getarg(1,ntile_out_str) endif -read(ntile_out_str,'(I2.2)') ntile_out +read(ntile_out_str,'(I1.1)') ntile_out open (unit=nlunit, file='input.nml', READONLY, status='OLD') read(nlunit,nam_stochy) close(nlunit) @@ -141,12 +139,10 @@ program standalone_stochy jec=Atm(1)%bd%jec nx=Atm(1)%npx-1 ny=Atm(1)%npy-1 -nx2=iec-isc+1 -ny2=jec-jsc+1 -allocate(workg(nx2,ny2)) -allocate(workg3d(nx2,ny2,nlevs)) -blksz=nx2 -nblks=ny2 +allocate(workg(nx,ny)) +allocate(workg3d(nx,ny,nlevs)) +nblks=ny +blksz=nx Allocate(Model%blksz(nblks)) Model%blksz(:)=blksz nthreads = omp_get_num_threads() @@ -184,26 +180,27 @@ program standalone_stochy Grid(j)%xlat(:)=Init_parm%xlat(:,j) Grid(j)%xlon(:)=Init_parm%xlon(:,j) enddo -allocate(grid_xt(nx2),grid_yt(ny2)) -do i=1,nx2 +allocate(grid_xt(nx),grid_yt(ny)) +do i=1,nx grid_xt(i)=i enddo -do i=1,ny2 +do i=1,ny grid_yt(i)=i enddo !setup GFS_coupling allocate(Coupling(nblks)) call init_stochastic_physics(Model, Init_parm, ntasks, nthreads) call get_outfile(fname) -write(strid,'(I2.1)') my_id+1 +write(strid,'(I1.1)') my_id+1 if (ntile_out.EQ.0) write_this_tile=.true. if ((my_id+1).EQ.ntile_out) write_this_tile=.true. -print*,trim(fname)//'.tile'//trim(adjustl(strid))//'.nc',write_this_tile +print*,trim(fname)//'.tile'//strid//'.nc',write_this_tile if (write_this_tile) then fid=30+my_id -ierr=nf90_create(trim(fname)//'.tile'//trim(adjustl(strid))//'.nc',cmode=NF90_CLOBBER,ncid=ncid) -ierr=NF90_DEF_DIM(ncid,"grid_xt",nx2,xt_dim_id) -ierr=NF90_DEF_DIM(ncid,"grid_yt",ny2,yt_dim_id) +!ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +ierr=NF90_DEF_DIM(ncid,"grid_xt",nx,xt_dim_id) +ierr=NF90_DEF_DIM(ncid,"grid_yt",ny,yt_dim_id) if (Model%do_skeb)ierr=NF90_DEF_DIM(ncid,"p_ref",nlevs,zt_dim_id) ierr=NF90_DEF_DIM(ncid,"time",NF90_UNLIMITED,time_dim_id) !> - Define the dimension variables. @@ -272,33 +269,33 @@ program standalone_stochy if (Model%do_skeb)allocate(Coupling(i)%skebu_wts(blksz,nlevs)) if (Model%do_skeb)allocate(Coupling(i)%skebv_wts(blksz,nlevs)) enddo -do i=1,10 +do i=1,200 Model%kdt=i ts=i/4.0 call run_stochastic_physics(Model, Grid, Coupling, nthreads) if (Model%me.EQ.0) print*,'sppt_wts=',i,Coupling(1)%sppt_wts(1,20) if (write_this_tile) then if (Model%do_sppt)then - do j=1,ny2 + do j=1,ny workg(:,j)=Coupling(j)%sppt_wts(:,20) enddo ierr=NF90_PUT_VAR(ncid,varid1,workg,(/1,1,i/)) endif if (Model%do_shum)then - do j=1,ny2 + do j=1,ny workg(:,j)=Coupling(j)%shum_wts(:,1) enddo ierr=NF90_PUT_VAR(ncid,varid2,workg,(/1,1,i/)) endif if (Model%do_skeb)then do k=1,nlevs - do j=1,ny2 + do j=1,ny workg3d(:,j,k)=Coupling(j)%skebu_wts(:,k) enddo enddo ierr=NF90_PUT_VAR(ncid,varid3,workg3d,(/1,1,1,i/)) do k=1,nlevs - do j=1,ny2 + do j=1,ny workg3d(:,j,k)=Coupling(j)%skebv_wts(:,k) enddo enddo diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index a09e5d9..9b89902 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -59,7 +59,6 @@ subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) nodes=ntasks gis_stochy%me=me gis_stochy%nodes=nodes -gis_stochy%yhalo=10 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 From 2d0836b5e3fd548e75a115583b00eb5d62ed0021 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 7 Apr 2020 19:52:37 +0000 Subject: [PATCH 08/10] added initialization of 0 for ca_deep array --- cellular_automata_global.F90 | 7 +++---- cellular_automata_sgs.F90 | 7 +++++++ 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index ff417f2..97f3224 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -88,6 +88,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & ! stop ! endif + call atmosphere_resolution (nlon, nlat, global=.false.) isize=nlon+2*halo jsize=nlat+2*halo @@ -129,7 +130,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & CA2(:,:) = 0.0 CA3(:,:) = 0.0 - !Put the blocks of model fields into a 2d array levs=nlev blocksz=blocksize @@ -169,7 +169,6 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo - !Initiate the cellular automaton with random numbers larger than nfracseed do j = 1,nyc @@ -188,8 +187,8 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & !we here set the "condition" variable to a different model field depending !on nf. (this is not used if ca_global = .true.) -do nf=1,nca !update each ca - + + do nf=1,nca !update each ca do j = 1,nyc do i = 1,nxc ilives_g(i,j)=int(real(nlives)*1.5*noise(i,j,nf)) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index 72b535f..ab4bbe0 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -522,6 +522,13 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & enddo enddo +if(kstep == 1)then +do j=1,nlat + do i=1,nlon + ca_plumes(i,j)=0. + enddo +enddo +endif !Put back into blocks 1D array to be passed to physics !or diagnostics output From 66b51ff88b617a56d3331a0461d8c8a0463f833f Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 29 Apr 2020 03:08:23 +0000 Subject: [PATCH 09/10] Updates to ensure reproducibility when option of seeding with new cell throughout the integration is activated --- cellular_automata_global.F90 | 2 +- cellular_automata_sgs.F90 | 2 +- update_ca.F90 | 47 ++++++++++++++++++++++++------------ 3 files changed, 34 insertions(+), 17 deletions(-) diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 index 97f3224..8ae6efc 100644 --- a/cellular_automata_global.F90 +++ b/cellular_automata_global.F90 @@ -203,7 +203,7 @@ subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & CA(:,:)=0. - call update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,iini_g,ilives_g, & + call update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,iseed_ca,CA,iini_g,ilives_g, & nlives, ncells, nfracseed, nseed,nthresh, nspinup,nf) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index ab4bbe0..96af296 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -363,7 +363,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & !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_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini,ilives, & + call update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,iseed_ca,CA,ca_plumes,iini,ilives, & nlives, ncells, nfracseed, nseed,nthresh,nspinup,nf,nca_plumes) if(nf==1)then diff --git a/update_ca.F90 b/update_ca.F90 index 77771ee..55d5305 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -16,12 +16,12 @@ module update_ca contains -subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,iini,ilives, & +subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,iseed_ca,CA,ca_plumes,iini,ilives, & nlives,ncells,nfracseed,nseed,nthresh,nspinup,nf,nca_plumes) implicit none -integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca +integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca,iseed_ca integer, intent(in) :: iini(nxc,nyc,nca) integer, intent(inout) :: ilives(nxc,nyc,nca) real, intent(out) :: CA(nlon,nlat) @@ -72,13 +72,20 @@ subroutine update_cells_sgs(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,ca_plumes,i 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 - !call mp_bcst(count5) - !call mp_bcst(count6) + if (iseed_ca == 0) then + ! generate a random seed from system clock and ens member number + call system_clock(count, count_rate, count_max) + ! iseed is elapsed time since unix epoch began (secs) + ! truncate to 4 byte integer + count_trunc = iscale*(count/iscale) + count5 = count - count_trunc + count6 = count5+9827 + else + ! don't rely on compiler to truncate integer(8) to integer(4) on + ! overflow, do wrap around explicitly. + count5 = mod(iseed_ca + 2147483648, 4294967296) - 2147483648 + count6 = count5 + 9827 + endif call random_setseed(count5) call random_number(noise1D1) @@ -298,12 +305,12 @@ end subroutine update_cells_sgs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,iini_g,ilives_g, & +subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,iseed_ca,CA,iini_g,ilives_g, & nlives,ncells,nfracseed,nseed,nthresh,nspinup,nf) implicit none -integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca +integer, intent(in) :: kstep,nxc,nyc,nlon,nlat,nxch,nych,nca,iseed_ca integer, intent(in) :: iini_g(nxc,nyc,nca), ilives_g(nxc,nyc) real, intent(out) :: CA(nlon,nlat) integer, intent(in) :: nlives, ncells, nseed, nspinup, nf @@ -350,10 +357,20 @@ subroutine update_cells_global(kstep,nca,nxc,nyc,nxch,nych,nlon,nlat,CA,iini_g,i 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 + if (iseed_ca == 0) then + ! generate a random seed from system clock and ens member number + call system_clock(count, count_rate, count_max) + ! iseed is elapsed time since unix epoch began (secs) + ! truncate to 4 byte integer + count_trunc = iscale*(count/iscale) + count5 = count - count_trunc + count6 = count5+9827 + else + ! don't rely on compiler to truncate integer(8) to integer(4) on + ! overflow, do wrap around explicitly. + count5 = mod(iseed_ca + 2147483648, 4294967296) - 2147483648 + count6 = count5 + 9827 + endif call random_setseed(count5) call random_number(noise1D1) From 05f564d7f9c2d58ab424e18aa61aec6ecc5601b3 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 1 May 2020 18:13:55 +0000 Subject: [PATCH 10/10] introducing k850 and k350 --- cellular_automata_sgs.F90 | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index 96af296..f859488 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -50,7 +50,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & integer :: inci, incj, nxc, nyc, nxch, nych integer :: halo, k_in, i, j, k, iec, jec, isc, jsc integer :: seed, ierr7,blk, ix, iix, count4,ih,jh -integer :: blocksz,levs +integer :: blocksz,levs,k350,k850 integer(8) :: count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 integer, allocatable :: iini(:,:,:),ilives(:,:,:),iini_g(:,:,:),ilives_g(:,:),ca_plumes(:,:) @@ -80,6 +80,18 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & halo=1 k_in=1 +if (nlev .EQ. 64) then + k350=29 + k850=13 +elseif (nlev .EQ. 127) then + k350=61 + k850=28 +else ! make a guess + k350=int(nlev/2) + k850=int(nlev/5) + print*,'this level selection is not supported, making an approximation for k350 and k850' +endif + nca_plumes = .true. !---------------------------------------------------------------------------- ! Get information about the compute domain, allocate fields on this @@ -180,8 +192,8 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & shalp(i,j) = Coupling(blk)%ca_shal(ix) gamt(i,j) = Coupling(blk)%ca_turb(ix) surfp(i,j) = Statein(blk)%pgr(ix) - humidity(i,j)=Statein(blk)%qgrs(ix,13,1) !about 850 hpa - do k = 1,29 !Lower troposphere: level 29 is about 350hPa + humidity(i,j)=Statein(blk)%qgrs(ix,k850,1) !about 850 hpa + do k = 1,k350 !Lower troposphere 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 @@ -198,11 +210,11 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i =1,nlon dp(i,j,1)=(surfp(i,j)-pressure(i,j,1)) - do k=2,29 + do k=2,k350 dp(i,j,k)=(pressure(i,j,k-1)-pressure(i,j,k)) enddo count1=0. - do k=1,29 + do k=1,k350 count1=count1+1. vertvelsum(i,j)=vertvelsum(i,j)+(omega(i,j,k)*dp(i,j,k)) enddo @@ -211,7 +223,7 @@ subroutine cellular_automata_sgs(kstep,Statein,Coupling,Diag,nblks,nlev, & do j=1,nlat do i=1,nlon - vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,29)) + vertvelmean(i,j)=vertvelsum(i,j)/(surfp(i,j)-pressure(i,j,k350)) enddo enddo