module cellular_automata_global_mod use update_ca, only : domain_global,iscnx_g,iecnx_g,jscnx_g,jecnx_g,isdnx_g,iednx_g,jsdnx_g,jednx_g, & nxncells_g,nyncells_g,csum,cold_start_ca_global implicit none contains subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cpl,ca3_cpl, & domain_in,nblks,isc,iec,jsc,jec,npx,npy,nlev, & nca,ncells,nlives,nfracseed,nseed,iseed_ca, mytile, & ca_smooth,nspinup,blocksize,nsmooth,ca_amplitude,mpiroot,mpicomm) use kinddef, only: kind_dbl_prec, kind_phys use update_ca, only: update_cells_global,define_ca_domain use halo_exchange, only: atmosphere_scalar_field_halo use random_numbers, only: random_01_CB use mpp_domains_mod, only: domain2D,mpp_get_global_domain,CENTER, mpp_get_data_domain, mpp_get_compute_domain,mpp_global_sum, & BITWISE_EFP_SUM, BITWISE_EXACT_SUM,mpp_define_io_domain,mpp_get_io_domain_layout use block_control_mod, only: block_control_type, define_blocks_packed use mpi_wrapper, only: mp_reduce_sum,mp_reduce_max,mp_reduce_min, & mpi_wrapper_initialize,mype,is_rootpe implicit none !L.Bengtsson, 2017-06 !P.Pegion, 2021-09 ! swtich to new random number generator and improve computational efficiency ! and remove unsued code. Also add restart capability ca_global !This program evolves a cellular automaton uniform over the globe integer, intent(in) :: kstep,ncells,nca,nlives,nseed,nspinup,nsmooth,mpiroot,mpicomm integer(kind=kind_dbl_prec), intent(in) :: iseed_ca integer, intent(in) :: mytile real(kind=kind_phys), intent(in) :: nfracseed,ca_amplitude logical, intent(in) :: ca_smooth,first_time_step, restart integer, intent(in) :: nblks,isc,iec,jsc,jec,npx,npy,nlev,blocksize real(kind=kind_phys), intent(out) :: ca1_cpl(:,:),ca2_cpl(:,:),ca3_cpl(:,:) type(domain2D), intent(inout) :: domain_in type(block_control_type) :: Atm_block integer :: nlon, nlat, isize,jsize,nf,nn integer :: inci, incj, nxc, nyc, nxch, nych integer :: halo, k_in, i, j, k integer :: seed, ierr7,blk, ix, iix, count4,ih,jh integer :: blocksz,levs integer,save :: isdnx,iednx,jsdnx,jednx integer,save :: iscnx,iecnx,jscnx,jecnx integer :: nxncells, nyncells integer(8) :: count, count_rate, count_max, count_trunc,nx_full integer(8) :: iscale = 10000000000 integer, allocatable :: iini_g(:,:,:),ilives_g(:,:) integer, allocatable :: io_layout(:) real(kind=kind_phys), allocatable :: field_out(:,:,:), field_smooth(:,:) real(kind=kind_phys), allocatable :: CA(:,:),CA1(:,:),CA2(:,:),CA3(:,:),CAprime(:,:) real*8 , allocatable :: noise(:,:,:) real*8 :: psum,CAmean,sq_diff,CAstdv,inv9 real*8 :: Detmax,Detmin logical,save :: block_message=.true. integer*8 :: i1,j1 integer :: ct !nca :: switch for number of cellular automata to be used. !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 if (nca .LT. 1) return ! Initialize MPI and OpenMP if (first_time_step) then call mpi_wrapper_initialize(mpiroot,mpicomm) end if halo=3 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 nlon=iec-isc+1 nlat=jec-jsc+1 isize=nlon+2*halo jsize=nlat+2*halo inci=ncells incj=ncells !--- get params from domain_ncellx for building board and board_halo !Get CA domain if(first_time_step)then if (.not. restart) then allocate(io_layout(2)) io_layout=mpp_get_io_domain_layout(domain_in) call define_ca_domain(domain_in,domain_global,halo,ncells,nxncells_g,nyncells_g) call mpp_define_io_domain(domain_global, io_layout) endif call mpp_get_data_domain (domain_global,isdnx_g,iednx_g,jsdnx_g,jednx_g) call mpp_get_compute_domain (domain_global,iscnx_g,iecnx_g,jscnx_g,jecnx_g) endif nxc = iecnx_g-iscnx_g+1 nyc = jecnx_g-jscnx_g+1 nxch = iednx_g-isdnx_g+1 nych = jednx_g-jsdnx_g+1 inv9=1.0/9.0 if(first_time_step) csum=int(6*(npx-1),kind=8)*int((npx-1),kind=8) !Allocate fields: allocate(field_out(isize,jsize,1)) allocate(field_smooth(nlon,nlat)) allocate(iini_g(nxc,nyc,nca)) allocate(ilives_g(nxc,nyc)) allocate(CA(nlon,nlat)) allocate(CAprime(nlon,nlat)) allocate(CA1(nlon,nlat)) allocate(CA2(nlon,nlat)) allocate(CA3(nlon,nlat)) allocate(noise(nxc,nyc,nca)) nx_full=int(npx-1,kind=8) !Initialize: noise(:,:,:) = 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 - can't use nlev and blocksize directly, !because the arguments to define_blocks_packed are intent(inout) and not intent(in). levs=nlev blocksz=blocksize call define_blocks_packed('cellular_automata', Atm_block, isc, iec, jsc, jec, levs, & blocksz, block_message) do j=1,nyc j1=j+(jsc-1)*ncells do i=1,nxc i1=i+(isc-1)*ncells 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 + mytile *( i1+nx_full*(j1-1)) ! no need to multply by 7 since time will be different in sgs else ! don't rely on compiler to truncate integer(8) to integer(4) on ! overflow, do wrap around explicitly. count4 = mod(((iseed_ca+7)*mytile)*(i1+nx_full*(j1-1))+ 2147483648, 4294967296) - 2147483648 endif ct=1 do nf=1,nca noise(i,j,nf)=real(random_01_CB(ct*kstep,count4),kind=8) ct=ct+1 enddo enddo enddo !Initiate the cellular automaton with random numbers larger than nfracseed do nf=1,nca 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,halo,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, & npx,npy,CA,iini_g,ilives_g, & nlives,ncells,nfracseed,nseed,nspinup,nf,mytile) if (ca_smooth) then field_out=0. field_out(1+halo:nlon+halo,1+halo:nlat+halo,1) = real(CA(1:nlon,1:nlat),kind=8) do nn=1,nsmooth !number of iterations for the smoothing. call atmosphere_scalar_field_halo(field_out,halo,isize,jsize,k_in,isc,iec,jsc,jec,npx,npy,domain_global) do j=1,nlat do i=1,nlon ih=i+halo jh=j+halo field_smooth(i,j)=(field_out(ih,jh,1)+field_out(ih-1,jh,1)+ & field_out(ih,jh-1,1)+field_out(ih+1,jh,1)+& field_out(ih,jh+1,1)+field_out(ih-1,jh-1,1)+& field_out(ih-1,jh+1,1)+field_out(ih+1,jh+1,1)+& field_out(ih+1,jh-1,1))*inv9 enddo enddo field_out(1+halo:nlon+halo,1+halo:nlat+halo,1) = field_smooth(1:nlon,1:nlat) enddo !nn do j=1,nlat do i=1,nlon CA(i,j)=field_smooth(i,j) enddo enddo endif !smooth !mean: !psum=SUM(CA) !call mp_reduce_sum(psum) psum= mpp_global_sum (domain_global, CA, flags=BITWISE_EXACT_SUM) CAmean=psum/csum !std: !sq_diff = 0. do j=1,nlat do i=1,nlon CAprime(i,j) = (CA(i,j)-CAmean)**2.0 enddo enddo !call mp_reduce_sum(sq_diff) sq_diff= mpp_global_sum (domain_global, CAprime, flags=BITWISE_EXACT_SUM) 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(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 ca1_cpl(blk,ix)=CA1(i,j) ca2_cpl(blk,ix)=CA2(i,j) ca3_cpl(blk,ix)=CA3(i,j) enddo enddo deallocate(field_out) deallocate(field_smooth) deallocate(iini_g) deallocate(ilives_g) deallocate(CA) deallocate(CAprime) deallocate(CA1) deallocate(CA2) deallocate(CA3) deallocate(noise) end subroutine cellular_automata_global end module cellular_automata_global_mod