From e2efcd4edf2b61883880ffbe2d7dbccfa7cef37c Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 13 Sep 2023 15:37:02 +0000 Subject: [PATCH 1/4] Add possibility to run the CA on a single tile, ensure decomposition of domain is proportional between the atmosphere domain and the higher resolution CA domain --- update_ca.F90 | 85 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 63 insertions(+), 22 deletions(-) diff --git a/update_ca.F90 b/update_ca.F90 index b55b5b1f..dc85e806 100644 --- a/update_ca.F90 +++ b/update_ca.F90 @@ -8,15 +8,15 @@ module update_ca use random_numbers, only: random_01_CB use mpi_wrapper, only: mype,mp_reduce_min,mp_reduce_max use mpp_domains_mod, only: domain2D,mpp_get_global_domain,CENTER, mpp_get_data_domain, mpp_get_compute_domain,mpp_get_ntile_count,& - mpp_define_mosaic,mpp_get_layout,mpp_define_io_domain,mpp_get_io_domain_layout + mpp_define_mosaic,mpp_get_layout,mpp_define_io_domain,mpp_get_io_domain_layout,mpp_get_domain_extents use mpp_mod, only: mpp_error, mpp_pe, mpp_root_pe, & NOTE, FATAL use fms2_io_mod, only: FmsNetcdfDomainFile_t, unlimited, & - open_file, close_file, & - register_axis, register_restart_field, & - register_variable_attribute, register_field, & - read_restart, write_restart, write_data, & - get_global_io_domain_indices, variable_exists + open_file, close_file, & + register_axis, register_restart_field, & + register_variable_attribute, register_field, & + read_restart, write_restart, write_data, & + get_global_io_domain_indices, variable_exists implicit none @@ -48,6 +48,7 @@ subroutine define_ca_domain(domain_in,domain_out,halo,ncells,nxncells_local,nync integer,intent(in) :: ncells,halo integer,intent(out) :: nxncells_local, nyncells_local integer :: layout(2) +integer, allocatable :: xextent(:,:), yextent(:,:) integer :: ntiles integer, allocatable :: pe_start(:), pe_end(:) @@ -59,7 +60,12 @@ subroutine define_ca_domain(domain_in,domain_out,halo,ncells,nxncells_local,nync call mpp_get_global_domain(domain_in,xsize=nx,ysize=ny,position=CENTER) call mpp_get_layout(domain_in,layout) ntiles = mpp_get_ntile_count(domain_in) - !write(1000+mpp_pe(),*) "nx,ny: ",nx,ny + allocate(xextent(layout(1),ntiles)) + allocate(yextent(layout(2),ntiles)) + call mpp_get_domain_extents(domain_in,xextent,yextent) + xextent=xextent*ncells + yextent=yextent*ncells + !write(1000+mpp_pe(),*) "nx,ny: ",nx,ny !write(1000+mpp_pe(),*) "layout: ",layout !--- define mosaic for domain_out refined by 'ncells' from domain_in @@ -72,9 +78,12 @@ subroutine define_ca_domain(domain_in,domain_out,halo,ncells,nxncells_local,nync pe_start(n) = mpp_root_pe() + (n-1)*layout(1)*layout(2) pe_end(n) = mpp_root_pe() + n*layout(1)*layout(2)-1 enddo - call define_cubic_mosaic(domain_out, nxncells_local-1, nyncells_local-1, layout, pe_start, pe_end, halo ) + call define_cubic_mosaic(domain_out, nxncells_local-1, nyncells_local-1, & + layout, pe_start, pe_end, halo, ntiles, xextent, yextent ) deallocate(pe_start) deallocate(pe_end) + deallocate(xextent) + deallocate(yextent) end subroutine define_ca_domain !--------------------------------------------------------------------------------------------- @@ -409,7 +418,7 @@ subroutine update_cells_sgs(kstep,halo,dt,initialize_ca,iseed_ca,first_flag,rest neighbours=0 birth=0 newcell=0 - + board_halo=0 !--- copy board into the halo-augmented board_halo board_halo(1+halo:nxc+halo,1+halo:nyc+halo,1) = real(board(1:nxc,1:nyc,1),kind=8) @@ -828,26 +837,29 @@ end subroutine update_cells_global ! and modified to make it simpler to use. ! domain_decomp in fv_mp_mod.F90 does something similar, but it does a ! few other unnecessary things (and requires more arguments). - subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo) + subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo, ntiles, xextent, yextent) type(domain2d), intent(inout) :: domain - integer, intent(in) :: ni, nj + integer, intent(in) :: ni, nj, ntiles, xextent(:,:), yextent(:,:) integer, intent(in) :: layout(:) integer, intent(in) :: pe_start(:), pe_end(:) integer, intent(in) :: halo !--- local variables - integer :: global_indices(4,6), layout2D(2,6) - integer, dimension(12) :: istart1, iend1, jstart1, jend1, tile1 - integer, dimension(12) :: istart2, iend2, jstart2, jend2, tile2 - integer :: ntiles, num_contact + integer,allocatable :: global_indices(:,:), layout2D(:,:) + integer,allocatable :: istart1(:), iend1(:), jstart1(:), jend1(:), tile1(:) + integer,allocatable :: istart2(:), iend2(:), jstart2(:), jend2(:), tile2(:) + integer :: num_contact integer :: i - ntiles = 6 - num_contact = 12 - if(size(pe_start(:)) .NE. 6 .OR. size(pe_end(:)) .NE. 6 ) call mpp_error(FATAL, & - "define_cubic_mosaic: size of pe_start and pe_end should be 6") if(size(layout) .NE. 2) call mpp_error(FATAL, & "define_cubic_mosaic: size of layout should be 2") - do i = 1, 6 + + if(ntiles==6)then + num_contact = 12 + allocate(global_indices(4,ntiles)) + allocate(layout2D(2,ntiles)) + allocate(istart1(num_contact), iend1(num_contact), jstart1(num_contact), jend1(num_contact), tile1(num_contact) ) + allocate(istart2(num_contact), iend2(num_contact), jstart2(num_contact), jend2(num_contact), tile2(num_contact) ) + do i = 1, ntiles layout2D(:,i) = layout(:) global_indices(1,i) = 1 global_indices(2,i) = ni @@ -914,11 +926,40 @@ subroutine define_cubic_mosaic(domain, ni, nj, layout, pe_start, pe_end, halo) istart1(12) = ni; iend1(12) = ni; jstart1(12) = 1; jend1(12) = nj istart2(12) = 1; iend2(12) = 1; jstart2(12) = 1; jend2(12) = nj + else if(ntiles==1) then !Single tile + + num_contact = 0 + + allocate(global_indices(4,ntiles)) + allocate(layout2D(2,ntiles)) + allocate(istart1(num_contact+1), iend1(num_contact+1), jstart1(num_contact+1), jend1(num_contact+1), tile1(num_contact+1) ) + allocate(istart2(num_contact+1), iend2(num_contact+1), jstart2(num_contact+1), jend2(num_contact+1), tile2(num_contact+1) ) + + do i = 1, ntiles + layout2D(:,i) = layout(:) + global_indices(1,i) = 1 + global_indices(2,i) = ni + global_indices(3,i) = 1 + global_indices(4,i) = nj + enddo + + else + + call mpp_error(FATAL, & + "ntiles should be either 6 or 1 to run cellular automata") + + endif !global or regional domain + call mpp_define_mosaic(global_indices, layout2D, domain, ntiles, & num_contact, tile1, tile2, istart1, iend1, jstart1, jend1, & istart2, iend2, jstart2, jend2, pe_start, pe_end, symmetry=.true., & - whalo=halo, ehalo=halo, shalo=halo, nhalo=halo, & - name='CA cubic mosaic') + whalo=halo, ehalo=halo, shalo=halo, nhalo=halo,name='CA cubic mosaic', & + xextent=xextent, yextent=yextent) + + deallocate(global_indices) + deallocate(layout2D) + deallocate(istart1, iend1, jstart1, jend1, tile1) + deallocate(istart2, iend2, jstart2, jend2, tile2) end subroutine define_cubic_mosaic From 4ba8fc1b5f9fb7bceecd0d069929b736c89463ad Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 13 Sep 2023 18:12:40 +0000 Subject: [PATCH 2/4] If running a regional model, don't mask out the CA outside of the tropical ocean --- cellular_automata_sgs.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/cellular_automata_sgs.F90 b/cellular_automata_sgs.F90 index d6aa7c2e..af2322d7 100644 --- a/cellular_automata_sgs.F90 +++ b/cellular_automata_sgs.F90 @@ -16,7 +16,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak use update_ca, only: update_cells_sgs, define_ca_domain 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_define_io_domain,mpp_get_io_domain_layout + mpp_define_io_domain,mpp_get_io_domain_layout,mpp_get_ntile_count use block_control_mod, only: block_control_type, define_blocks_packed use time_manager_mod, only: time_type use mpi_wrapper, only: mype,mp_reduce_max, & @@ -75,7 +75,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak logical :: nca_plumes logical,save :: first_flag integer*8 :: i1,j1 -integer :: ct +integer :: ct,ntiles real :: dz,invgrav !nca :: switch for number of cellular automata to be used. @@ -124,7 +124,8 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak !Set time and length scales: call mpp_get_global_domain(domain_in,xsize=nx,ysize=ny,position=CENTER) - + ntiles = mpp_get_ntile_count(domain_in) + if(mype == 1)then write(*,*)'ncells=',ncells write(*,*)'nlives=',nlives @@ -348,6 +349,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak !Limit CA activity to the Tropical Ocean +if(ntiles==6)then do j=1,nlat do i=1,nlon if(ssti(i,j) < 300. .or. lsmski(i,j) /= 0. .or. lakei(i,j) > 0.0)then @@ -355,6 +357,7 @@ subroutine cellular_automata_sgs(kstep,dtf,restart,first_time_step,sst,lsmsk,lak endif enddo enddo +endif !Put back into blocks 1D array to be passed to physics !or diagnostics output From 6905288dbedcddad1fe3c5f2b0a4dd6543421828 Mon Sep 17 00:00:00 2001 From: Bing Fu <48262811+bingfu-NOAA@users.noreply.github.com> Date: Fri, 15 Sep 2023 18:01:17 +0000 Subject: [PATCH 3/4] Adding pbl_taper to namelist to define PBL tapering profile (#65) Co-authored-by: bing fu --- compns_stochy.F90 | 5 +++-- stochastic_physics.F90 | 5 +++-- stochy_namelist_def.F90 | 1 + 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index c18ea839..3f3704ab 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -62,7 +62,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper namelist /nam_sfcperts/lndp_type,lndp_model_type, lndp_var_list, lndp_prt_list, & iseed_lndp, lndp_tau,lndp_lscale ! For SPP physics parameterization perterbations @@ -144,6 +144,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) spp_sigtop2 = 0.025 ! reduce amplitude of sppt near surface (lowest 2 levels) sppt_sfclimit = .false. + pbl_taper = (/0.0,0.5,1.0,1.0,1.0,1.0,1.0/) ! gaussian or power law variance spectrum for skeb (0: gaussian, 1: ! power law). If power law, skeb_lscale interpreted as a power not a ! length scale. @@ -444,7 +445,7 @@ subroutine compns_stochy_ocn (deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper namelist /nam_sfcperts/lndp_type,lndp_model_type,lndp_var_list, lndp_prt_list, iseed_lndp, & lndp_tau,lndp_lscale diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 18ff0a2d..506f8fb6 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -162,8 +162,9 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in endif enddo if (sppt_sfclimit) then - vfact_sppt(2)=vfact_sppt(3)*0.5 - vfact_sppt(1)=0.0 + do k=1,7 + vfact_sppt(k)=pbl_taper(k) + enddo endif if (is_rootpe()) then do k=1,levs diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index bb8f1eda..08b9dfd9 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -27,6 +27,7 @@ module stochy_namelist_def real(kind=kind_phys), dimension(5) :: shum,shum_lscale,shum_tau real(kind=kind_phys), dimension(5) :: epbl,epbl_lscale,epbl_tau real(kind=kind_phys), dimension(5) :: ocnsppt,ocnsppt_lscale,ocnsppt_tau + real(kind=kind_phys), dimension(7) :: pbl_taper integer,dimension(5) ::skeb_vfilt integer(kind=8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_epbl2 logical stochini,sppt_logit,new_lscale From 2dcefcecef21f465ee0c932b467372a63dd451f6 Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Fri, 15 Sep 2023 13:22:14 -0600 Subject: [PATCH 4/4] Revert "Adding pbl_taper to namelist to define PBL tapering profile (#65)" This reverts commit 6905288dbedcddad1fe3c5f2b0a4dd6543421828. --- compns_stochy.F90 | 5 ++--- stochastic_physics.F90 | 5 ++--- stochy_namelist_def.F90 | 1 - 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index 3f3704ab..c18ea839 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -62,7 +62,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt namelist /nam_sfcperts/lndp_type,lndp_model_type, lndp_var_list, lndp_prt_list, & iseed_lndp, lndp_tau,lndp_lscale ! For SPP physics parameterization perterbations @@ -144,7 +144,6 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) spp_sigtop2 = 0.025 ! reduce amplitude of sppt near surface (lowest 2 levels) sppt_sfclimit = .false. - pbl_taper = (/0.0,0.5,1.0,1.0,1.0,1.0,1.0/) ! gaussian or power law variance spectrum for skeb (0: gaussian, 1: ! power law). If power law, skeb_lscale interpreted as a power not a ! length scale. @@ -445,7 +444,7 @@ subroutine compns_stochy_ocn (deltim,iret) skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale, & epbl,epbl_lscale,epbl_tau,iseed_epbl, & - ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt,pbl_taper + ocnsppt,ocnsppt_lscale,ocnsppt_tau,iseed_ocnsppt namelist /nam_sfcperts/lndp_type,lndp_model_type,lndp_var_list, lndp_prt_list, iseed_lndp, & lndp_tau,lndp_lscale diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 506f8fb6..18ff0a2d 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -162,9 +162,8 @@ subroutine init_stochastic_physics(levs, blksz, dtp, sppt_amp, input_nml_file_in endif enddo if (sppt_sfclimit) then - do k=1,7 - vfact_sppt(k)=pbl_taper(k) - enddo + vfact_sppt(2)=vfact_sppt(3)*0.5 + vfact_sppt(1)=0.0 endif if (is_rootpe()) then do k=1,levs diff --git a/stochy_namelist_def.F90 b/stochy_namelist_def.F90 index 08b9dfd9..bb8f1eda 100644 --- a/stochy_namelist_def.F90 +++ b/stochy_namelist_def.F90 @@ -27,7 +27,6 @@ module stochy_namelist_def real(kind=kind_phys), dimension(5) :: shum,shum_lscale,shum_tau real(kind=kind_phys), dimension(5) :: epbl,epbl_lscale,epbl_tau real(kind=kind_phys), dimension(5) :: ocnsppt,ocnsppt_lscale,ocnsppt_tau - real(kind=kind_phys), dimension(7) :: pbl_taper integer,dimension(5) ::skeb_vfilt integer(kind=8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb,iseed_epbl,iseed_ocnsppt,iseed_epbl2 logical stochini,sppt_logit,new_lscale