diff --git a/cellular_automata_global.F90 b/cellular_automata_global.F90 new file mode 100644 index 0000000..8ae6efc --- /dev/null +++ b/cellular_automata_global.F90 @@ -0,0 +1,346 @@ +subroutine cellular_automata_global(kstep,Statein,Coupling,Diag,nblks,nlev, & + nca,ncells,nlives,nfracseed,nseed,nthresh,ca_global,ca_sgs,iseed_ca, & + ca_smooth,nspinup,blocksize,nsmooth,ca_amplitude) + +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,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) +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_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 +real(kind=kind_phys) :: Detmax(nca),Detmin(nca),Detmean(nca) +logical,save :: block_message=.true. + + +!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 + + +halo=1 +k_in=1 + +!---------------------------------------------------------------------------- +! Get information about the compute domain, allocate fields on this +! domain + +! WRITE(*,*)'Entering cellular automata calculations' + +! Some security checks for namelist combinations: + if(nca > 3)then + write(0,*)'Namelist option nca cannot be larger than 3 - 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(field_in(nlon*nlat,1)) + allocate(field_out(isize,jsize,1)) + allocate(field_smooth(nlon,nlat)) + allocate(iini_g(nxc,nyc,nca)) + allocate(ilives_g(nxc,nyc)) + allocate(vertvelhigh(nxc,nyc)) + allocate(CA(nlon,nlat)) + allocate(CA1(nlon,nlat)) + allocate(CA2(nlon,nlat)) + allocate(CA3(nlon,nlat)) + allocate(noise(nxc,nyc,nca)) + allocate(noise1D(nxc*nyc)) + + !Initialize: + + noise(:,:,:) = 0.0 + noise1D(:) = 0.0 + iini_g(:,:,:) = 0 + ilives_g(:,:) = 0 + CA1(:,:) = 0.0 + CA2(:,:) = 0.0 + CA3(:,:) = 0.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 + +!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,iseed_ca,CA,iini_g,ilives_g, & + nlives, ncells, nfracseed, nseed,nthresh, nspinup,nf) + + +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)=(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))/18. + enddo +enddo + +do j=1,nlat + do i=1,nlon + CA(i,j)=field_smooth(i,j) + enddo +enddo + +enddo !nn +endif !smooth + +!!!!Post processing, should be made into a separate subroutine + +Detmax(1)=maxval(CA) +call mp_reduce_max(Detmax(1)) +Detmin(1)=minval(CA) +call mp_reduce_min(Detmin(1)) + +do j=1,nlat + do i=1,nlon + CA(i,j) = ((CA(i,j) - Detmin(1))/(Detmax(1)-Detmin(1))) + enddo +enddo + +!mean: +CAmean=0. +psum=0. +csum=0. +do j=1,nlat + do i=1,nlon + psum=psum+(CA(i,j)) + csum=csum+1 + enddo +enddo + +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) + +CAmean=psum/csum + +!std: +CAstdv=0. +sq_diff = 0. +do j=1,nlat + do i=1,nlon + sq_diff = sq_diff + (CA(i,j)-CAmean)**2.0 + enddo +enddo + +call mp_reduce_sum(sq_diff) + +CAstdv = sqrt(sq_diff/csum) + +!Transform to mean of 1 and ca_amplitude standard deviation + +do j=1,nlat + do i=1,nlon + CA(i,j)=1.0 + (CA(i,j)-CAmean)*(ca_amplitude/CAstdv) + enddo +enddo + +do j=1,nlat + do i=1,nlon + 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 +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) + 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(field_in) + deallocate(field_out) + deallocate(field_smooth) + deallocate(iini_g) + deallocate(ilives_g) + deallocate(CA) + 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 53% rename from cellular_automata.f90 rename to cellular_automata_sgs.F90 index b8da0c2..f859488 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,93 +46,85 @@ 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,k350,k850 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(:,:),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,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,condmax 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 -nca_plumes = .false. +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 ! 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: + levs=nlev allocate(cloud(nlon,nlat)) - allocate(omega(nlon,nlat,k350)) - allocate(pressure(nlon,nlat,k350)) + allocate(omega(nlon,nlat,levs)) + allocate(pressure(nlon,nlat,levs)) allocate(humidity(nlon,nlat)) - allocate(dp(nlon,nlat,k350)) + allocate(uwind(nlon,nlat)) + allocate(dp(nlon,nlat,levs)) allocate(rho(nlon,nlat)) allocate(surfp(nlon,nlat)) allocate(vertvelmean(nlon,nlat)) @@ -137,44 +133,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(conditiongrid(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. + conditiongrid(:,:)=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,22 +181,29 @@ 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) + 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) humidity(i,j)=Statein(blk)%qgrs(ix,k850,1) !about 850 hpa - do k = 1,k350 !Lower troposphere: level k350 is about 350hPa + do k = 1,k350 !Lower troposphere 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. @@ -235,26 +236,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 +264,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)=conditiongrid(inci/ncells,incj/ncells) !conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -293,9 +291,12 @@ subroutine cellular_automata(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)=int(real(nlives)*alpha*noise(i,j,nf)) + do i = 1,nxc + ilives(i,j,nf)=real(nlives)*(condition(i,j)/condmax) enddo enddo @@ -304,7 +305,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)=conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -315,40 +316,22 @@ subroutine cellular_automata(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)=int(real(nlives)*alpha*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)=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 + else - do j = 1,nyc - do i = 1,nxc - ilives(i,j)=int(real(nlives)*alpha*noise(i,j,nf)) - enddo - enddo - - elseif(nf==4)then inci=ncells incj=ncells do j=1,nyc do i=1,nxc - condition(i,j)=cape(inci/ncells,incj/ncells) + condition(i,j)=conditiongrid(inci/ncells,incj/ncells) if(i.eq.inci)then inci=inci+ncells endif @@ -359,38 +342,18 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif enddo - do j = 1,nyc - do i = 1,nxc - ilives(i,j)=int(real(nlives)*alpha*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)=int(real(nlives)*alpha*noise(i,j,nf)) + ilives(i,j,nf)=real(nlives)*(condition(i,j)/condmax) enddo enddo - endif !nf - + !Vertical velocity has its own variable in order to condition on combination !of "condition" and vertical velocity. @@ -409,144 +372,139 @@ subroutine cellular_automata(kstep,Statein,Coupling,Diag,nblks,nlev, & endif 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 + + 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) - 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, & - 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 -do j=1,nlat - do i=1,nlon - field_in(i+(j-1)*nlon,1)=CAavg(i,j) - 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. -!do j=1,nlat -! do i=1,nlon -! psum=psum+CAavg(i,j) -! csum=csum+1 -! enddo -!enddo +Detmin(1) = minval(CA_DEEP,CA_DEEP.NE.0) +call mp_reduce_min(Detmin(1)) -!call mp_reduce_sum(psum) -!call mp_reduce_sum(csum) +!Shallow convection ============================================================ -!CAmean=psum/csum +!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)) -!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 +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 -!call mp_reduce_sum(sq_diff) +!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 -!CAstdv = sqrt( sq_diff / (csum-1.0) ) +call mp_reduce_sum(psum) +call mp_reduce_sum(csum) -!do j=1,nlat -! do i=1,nlon -! CAavg(i,j)=(CAavg(i,j)-CAmean)/CAstdv -! enddo -!enddo +CAmean=psum/csum -!endif +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 the range for the nca individual ca_sgs patterns: -if(ca_sgs)then +!Turbulence ============================================================================= -Detmax(1)=maxval(CA_DEEP(:,:)) -call mp_reduce_max(Detmax(1)) +!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_DEEP(i,j)>0.)then - CA_DEEP(i,j)=CA_DEEP(i,j)/Detmax(1) !Now the range goes from 0-1 - endif + 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,54 +517,55 @@ 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 -!!! - +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 + if(conditiongrid(i,j) == 0)then + CA_DEEP(i,j)=0. + ca_plumes(i,j)=0. + endif 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 - + 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_deep(ix)=ca_plumes(i,j) + Diag(blk)%ca_turb(ix)=conditiongrid(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_deep(ix)=ca_plumes(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(conditiongrid) + deallocate(shalp) + deallocate(gamt) deallocate(rho) deallocate(surfp) deallocate(vertvelmean) @@ -620,14 +579,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 6a7575c..8097503 100644 --- a/makefile +++ b/makefile @@ -22,8 +22,6 @@ FFLAGS += -I../FV3/gfsphysics/ -I../FV3/atmos_cubed_sphere -I$(FMS_DIR) -I../F 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/plumes.f90 b/plumes.f90 index 410cc47..80b301f 100644 --- a/plumes.f90 +++ b/plumes.f90 @@ -1,31 +1,24 @@ -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" - -!! 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 -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 new file mode 100644 index 0000000..55d5305 --- /dev/null +++ b/update_ca.F90 @@ -0,0 +1,533 @@ +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,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,iseed_ca +integer, intent(in) :: iini(nxc,nyc,nca) +integer, intent(inout) :: ilives(nxc,nyc,nca) +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(:),B(:) +integer,allocatable :: AG(:,:) +integer :: inci, incj, i, j, k, iii,sub,spinup,it,halo,k_in,isize,jsize +integer :: ih, jh,kend +real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh +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 + +integer(8) :: count, count_rate, count_max, count_trunc +integer(8) :: iscale = 10000000000 +integer :: count5, count6 +type(random_stat) :: rstate +real, dimension(nxc,nyc) :: NOISE_A, NOISE_B +real, dimension(nxc*nyc) :: noise1D2, noise1D1 + + +!------------------------------------------------------------------------------------------------- +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(field_in))then + allocate(field_in(nxc*nyc,1)) + endif + if(.not. allocated(board_halo))then + allocate(board_halo(nxch,nych,1)) + endif + + noise1D1 = 0.0 + noise1D2 = 0.0 + + 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) + + call random_setseed(count6) + call random_number(noise1D2) + + 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 + + newseed = 0 + if(mod(kstep,nseed) == 0 .and. kstep >= 2)then + do j=1,nyc + do i=1,nxc + 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 + + 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 + field_in=0 + maxlives = 0 + maxliveshigh =0 + + +!Step 4 - Compute the neighbourhood + + 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 + +! Step 5 - Check rules; + + !birth + + do j=1,nyc + do i=1,nxc + 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 + 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) + 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 + if (.not. allocated(B))then + allocate(B(kend)) + endif + if (.not. allocated(AG))then + allocate(AG(ncells,ncells)) + endif + + ca_plumes(:,:)=0 + inci=ncells + incj=ncells + 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,AG,onegrid,ncells,ncells,kend) + do k=1,kend + 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 + 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,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,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 +real, intent(in) :: nfracseed, nthresh +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 :: ih, jh +real :: threshc,threshk,wp_max,wp_min,mthresh,kthresh +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(8) :: count, count_rate, count_max, count_trunc +integer(8) :: iscale = 10000000000 +integer :: count5, count6 +type(random_stat) :: rstate +real, dimension(nxc,nyc) :: NOISE_A, NOISE_B, g2D +real, dimension(nxc*nyc) :: noise1D2, noise1D1 + + +!------------------------------------------------------------------------------------------------- +halo=1 +isize=nlon+2*halo +jsize=nlat+2*halo +k_in=1 + + 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(field_in))then + allocate(field_in(nxc*nyc,1)) + endif + if(.not. allocated(board_halo))then + allocate(board_halo(nxch,nych,1)) + endif + + + !random numbers: + noise1D1 = 0.0 + noise1D2 = 0.0 + + 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) + + 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 + 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 + +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