Skip to content

2D advection of cellular automata #62

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pjpegion merged 4 commits into from
Sep 5, 2023
Merged

2D advection of cellular automata #62

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
696b778
2D advection of cellular automata
lisa-bengtsson Aug 16, 2023
936fdb1
Merge branch 'master' of https://github.com/NOAA-PSL/stochastic_physi…
lisa-bengtsson Aug 23, 2023
f6994cc
Address review comment - invgrav
lisa-bengtsson Aug 23, 2023
8015b8f
merge upstream
lisa-bengtsson Aug 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 10 additions & 5 deletions cellular_automata_global.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp
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
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
Expand Down Expand Up @@ -50,6 +50,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp
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(:,:,:)
Expand All @@ -73,7 +74,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp
call mpi_wrapper_initialize(mpiroot,mpicomm)
end if

halo=1
halo=3
k_in=1

!----------------------------------------------------------------------------
Expand Down Expand Up @@ -102,8 +103,12 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp
!Get CA domain

if(first_time_step)then
! if (.not. restart) call define_ca_domain(domain_in,domain_global,ncells,nxncells_g,nyncells_g)
domain_global=domain_in
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
Expand Down Expand Up @@ -203,7 +208,7 @@ subroutine cellular_automata_global(kstep,restart,first_time_step,ca1_cpl,ca2_cp

CA(:,:)=0.

call update_cells_global(kstep,first_time_step,iseed_ca,restart,nca,nxc,nyc,nxch,nych,nlon,nlat,isc,iec,jsc,jec, &
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)

Expand Down
86 changes: 68 additions & 18 deletions cellular_automata_sgs.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@ module cellular_automata_sgs_mod

contains

subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lake,condition_cpl, &
subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lake,uwind,vwind,height,dx,condition_cpl, &
ca_deep_cpl,ca_turb_cpl,ca_shal_cpl,domain_in, &
nblks,isc,iec,jsc,jec,npx,npy,nlev,nthresh, mytile, &
nca,ncells,nlives,nfracseed,nseed,iseed_ca, &
nca,ncells,nlives,nfracseed,nseed,iseed_ca,ca_advect, &
nspinup,ca_trigger,blocksize,mpiroot,mpicomm)

use kinddef, only: kind_phys,kind_dbl_prec
Expand Down Expand Up @@ -37,16 +37,19 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
! swtich to new random number generator and improve computational efficiency
! and remove unsued code

!L.Bengtsson, 2023-05
!Add horizontal advection of CA cells

!This routine produces an output field CA_DEEP for coupling to convection (saSAS).
!CA_DEEP can be either number of plumes in a cluster (nca_plumes=true) or updraft
!area fraction (nca_plumes=false)

integer,intent(in) :: kstep,ncells,nca,nlives,nseed,nspinup,mpiroot,mpicomm,mytile
integer(kind=kind_dbl_prec), intent(in) :: iseed_ca
real(kind=kind_phys), intent(in) :: nfracseed,dtf,nthresh
logical,intent(in) :: restart,ca_trigger,first_time_step
logical,intent(in) :: restart,ca_trigger,first_time_step,ca_advect
integer, intent(in) :: nblks,isc,iec,jsc,jec,npx,npy,nlev,blocksize
real(kind=kind_phys), intent(in) :: sst(:,:),lsmsk(:,:),lake(:,:)
real(kind=kind_phys), intent(in) :: sst(:,:),lsmsk(:,:),lake(:,:),uwind(:,:,:),dx(:,:),vwind(:,:,:),height(:,:,:)
real(kind=kind_phys), intent(inout) :: condition_cpl(:,:)
real(kind=kind_phys), intent(inout) :: ca_deep_cpl(:,:)
real(kind=kind_phys), intent(inout) :: ca_turb_cpl(:,:)
Expand All @@ -58,21 +61,22 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
integer :: inci, incj, nxc, nyc, nxch, nych, nx, ny
integer :: halo, k_in, i, j, k
integer :: seed, ierr7,blk, ix, iix, count4,ih,jh
integer :: blocksz,levs
integer :: blocksz,levs,u200,u850
integer, save :: initialize_ca
integer(8) :: count, count_rate, count_max, count_trunc,nx_full
integer(8) :: iscale = 10000000000
integer, allocatable :: iini(:,:,:),ilives_in(:,:,:),ca_plumes(:,:),io_layout(:)
real(kind=kind_phys), allocatable :: ssti(:,:),lsmski(:,:),lakei(:,:)
real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),conditiongrid(:,:)
real(kind=kind_phys), allocatable :: CA_DEEP(:,:)
real(kind=kind_phys), allocatable :: ssti(:,:),lsmski(:,:),lakei(:,:),uwindi(:,:),vwindi(:,:),dxi(:,:),heighti(:,:,:)
real(kind=kind_phys), allocatable :: CA(:,:),condition(:,:),uhigh(:,:),vhigh(:,:),dxhigh(:,:),conditiongrid(:,:)
real(kind=kind_phys), allocatable :: CA_DEEP(:,:),zi(:,:,:),sumx(:,:)
real*8 , allocatable :: noise(:,:,:)
real(kind=kind_phys) :: condmax,condmaxinv,livesmax,livesmaxinv,factor,dx,pi,re
real(kind=kind_phys) :: condmax,condmaxinv,livesmax,livesmaxinv,factor,pi,re
logical,save :: block_message=.true.
logical :: nca_plumes
logical,save :: first_flag
integer*8 :: i1,j1
integer :: ct
real :: dz,grav

!nca :: switch for number of cellular automata to be used.
! :: for the moment only 1 CA can be used
Expand All @@ -88,8 +92,13 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
call mpi_wrapper_initialize(mpiroot,mpicomm)
end if

halo=1
halo=3
k_in=1
!Right now these values are experimental:
u200=56
u850=13
!Gravitational acceleration:
grav=9.81

nca_plumes = .true.

Expand Down Expand Up @@ -132,15 +141,12 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
if (.not.restart) then
allocate(io_layout(2))
io_layout=mpp_get_io_domain_layout(domain_in)
call define_ca_domain(domain_in,domain_sgs,ncells,nxncells,nyncells)
call define_ca_domain(domain_in,domain_sgs,halo,ncells,nxncells,nyncells)
call mpp_define_io_domain(domain_sgs, io_layout)
endif
call mpp_get_data_domain (domain_sgs,isdnx,iednx,jsdnx,jednx)
call mpp_get_compute_domain (domain_sgs,iscnx,iecnx,jscnx,jecnx)
!write(1000+mpp_pe(),*) "nxncells,nyncells: ",nxncells,nyncells
!write(1000+mpp_pe(),*) "iscnx,iecnx,jscnx,jecnx: ",iscnx,iecnx,jscnx,jecnx
!write(1000+mpp_pe(),*) "isdnx,iednx,jsdnx,jednx: ",isdnx,iednx,jsdnx,jednx
endif
endif

nxc = iecnx-iscnx+1
nyc = jecnx-jscnx+1
Expand All @@ -152,9 +158,18 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
allocate(ssti(nlon,nlat))
allocate(lsmski(nlon,nlat))
allocate(lakei(nlon,nlat))
allocate(uwindi(nlon,nlat))
allocate(vwindi(nlon,nlat))
allocate(heighti(nlon,nlat,nlev))
allocate(zi(nlon,nlat,nlev))
allocate(sumx(nlon,nlat))
allocate(dxi(nlon,nlat))
allocate(iini(nxc,nyc,nca))
allocate(ilives_in(nxc,nyc,nca))
allocate(condition(nxc,nyc))
allocate(uhigh(nxc,nyc))
allocate(dxhigh(nxc,nyc))
allocate(vhigh(nxc,nyc))
allocate(conditiongrid(nlon,nlat))
allocate(CA(nlon,nlat))
allocate(ca_plumes(nlon,nlat))
Expand All @@ -163,6 +178,9 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
!Initialize:
ilives_in(:,:,:) = 0
iini(:,:,:) = 0
sumx(:,:) = 0.
uwindi(:,:) = 0.
vwindi(:,:) =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).
Expand All @@ -172,6 +190,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
call define_blocks_packed('cellular_automata', Atm_block, isc, iec, jsc, jec, levs, &
blocksz, block_message)


do blk = 1,Atm_block%nblks
do ix = 1, Atm_block%blksz(blk)
i = Atm_block%index(blk)%ii(ix) - isc + 1
Expand All @@ -180,13 +199,39 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
ssti(i,j) = sst(blk,ix)
lsmski(i,j) = lsmsk(blk,ix)
lakei(i,j) = lake(blk,ix)
dxi(i,j) = dx(blk,ix)
do k = 1,levs
heighti(i,j,k) = height(blk,ix,k)/grav
Comment thread
lisa-bengtsson marked this conversation as resolved.
Outdated
enddo
do k = 2,levs
zi(i,j,k)=0.5*(heighti(i,j,k)+heighti(i,j,k-1))
enddo
do k = u850,u200
dz=zi(i,j,k)-zi(i,j,k-1)
uwindi(i,j) = uwindi(i,j) + uwind(blk,ix,k)*dz
vwindi(i,j) = vwindi(i,j) + vwind(blk,ix,k)*dz
sumx(i,j) = sumx(i,j) + dz
enddo
enddo
enddo

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
uwindi(i,j)=uwindi(i,j)/sumx(i,j)
vwindi(i,j)=vwindi(i,j)/sumx(i,j)
enddo
enddo


!Initialize the CA when the condition field is populated
do j=1,nyc
do i=1,nxc
condition(i,j)=conditiongrid(inci/ncells,incj/ncells)
uhigh(i,j)=uwindi(inci/ncells,incj/ncells)
vhigh(i,j)=vwindi(inci/ncells,incj/ncells)
dxhigh(i,j)=dxi(inci/ncells,incj/ncells)/ncells !dx on the finer grid
if(i.eq.inci)then
inci=inci+ncells
endif
Expand Down Expand Up @@ -277,9 +322,9 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak

!Calculate neighbours and update the automata
do nf=1,nca
call update_cells_sgs(kstep,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc, &
nxch,nych,nlon,nlat,isc,iec,jsc,jec, &
npx,npy,CA,ca_plumes,iini,ilives_in, &
call update_cells_sgs(kstep,halo,dtf,initialize_ca,iseed_ca,first_flag,restart,first_time_step,nca,nxc,nyc, &
nxch,nych,nlon,nlat,isc,iec,jsc,jec,ca_advect, &
npx,npy,CA,ca_plumes,iini,ilives_in,uhigh,vhigh,dxhigh, &
nlives,nfracseed,nseed,nspinup,nf,nca_plumes,ncells,mytile)

if(nca_plumes)then
Expand Down Expand Up @@ -328,6 +373,11 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak
deallocate(lakei)
deallocate(iini)
deallocate(ilives_in)
deallocate(uwindi)
deallocate(vwindi)
deallocate(uhigh)
deallocate(vhigh)
deallocate(dxhigh)
deallocate(condition)
deallocate(CA)
deallocate(ca_plumes)
Expand Down
4 changes: 2 additions & 2 deletions halo_exchange.fv3.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,8 @@ subroutine atmosphere_scalar_field_halo (data, halo, isize, jsize, ksize, &
if (halo == 1) then
call mpp_update_domains(data, domain_for_coupler, flags=mpp_flags, complete=.true.)
! Not needed for cellular automata code
!elseif (halo == 3) then
! call mpp_update_domains(data, Atm(mytile)%domain, flags=mpp_flags, complete=.true.)
elseif (halo == 3) then
call mpp_update_domains(data, domain_for_coupler, flags=mpp_flags, complete=.true.)
else
call mpp_error(FATAL, modname//' - unsupported halo size')
endif
Expand Down
Loading