Skip to content

add pbl tapering for SPPT #68

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 6 commits into from
Sep 26, 2023
Merged

add pbl tapering for SPPT #68

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
97388fd
Adding pbl_taper to namelist to define PBL tapering profile
Sep 9, 2023
e2efcd4
Add possibility to run the CA on a single tile, ensure decomposition …
lisa-bengtsson Sep 13, 2023
4ba8fc1
If running a regional model, don't mask out the CA outside of the tro…
lisa-bengtsson Sep 13, 2023
6d57dfc
Merge branch 'NOAA-PSL:master' into master
bingfu-NOAA Sep 15, 2023
8c4433e
Merge branch 'master' of https://github.com/NOAA-PSL/stochastic_physi…
lisa-bengtsson Sep 18, 2023
db232e0
Merge pull request #1 from lisa-bengtsson/CA_in_rrfs
bingfu-NOAA Sep 26, 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
9 changes: 6 additions & 3 deletions cellular_automata_sgs.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_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, &
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -348,13 +349,15 @@ 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
CA_DEEP(i,j)=0.
endif
enddo
enddo
endif

!Put back into blocks 1D array to be passed to physics
!or diagnostics output
Expand Down
5 changes: 3 additions & 2 deletions compns_stochy.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down
5 changes: 3 additions & 2 deletions stochastic_physics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions stochy_namelist_def.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
85 changes: 63 additions & 22 deletions update_ca.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(:)

Expand All @@ -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
Expand All @@ -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
!---------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down