Skip to content

Commit

Permalink
Ca develop (NOAA-EMC#87)
Browse files Browse the repository at this point in the history 
Updates to cellular automata for deep convection. Some cleaning. Changes goes with updates in stochastic_physics and ccpp/physics repositories.
  • Loading branch information
@lisa-bengtsson
lisa-bengtsson authored May 4, 2020
1 parent bc9aa3b commit e8fd6ee
Show file tree
Hide file tree
Showing 10 changed files with 353 additions and 213 deletions.
41 changes: 27 additions & 14 deletions atmos_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,14 +312,18 @@ subroutine update_atmos_radiation_physics (Atmos)
end if

if(IPD_Control%do_ca)then
! DH* The current implementation of cellular_automata assumes that all blocksizes are the
! same, this is tested in the initialization call to cellular_automata, no need to redo *DH
call cellular_automata(IPD_Control%kdt, IPD_Data(:)%Statein, IPD_Data(:)%Coupling, IPD_Data(:)%Intdiag, &
Atm_block%nblks, IPD_Control%levs, IPD_Control%nca, IPD_Control%ncells, &
IPD_Control%nlives, IPD_Control%nfracseed, IPD_Control%nseed, &
IPD_Control%nthresh, IPD_Control%ca_global, IPD_Control%ca_sgs, &
IPD_Control%iseed_ca, IPD_Control%ca_smooth, IPD_Control%nspinup, &
Atm_block%blksz(1))
if(IPD_Control%ca_sgs == .true.)then
call cellular_automata_sgs(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
IPD_Control%nca,IPD_Control%ncells,IPD_Control%nlives,IPD_Control%nfracseed,&
IPD_Control%nseed,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1))
endif
if(IPD_Control%ca_global == .true.)then
call cellular_automata_global(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
IPD_Control%nca_g,IPD_Control%ncells_g,IPD_Control%nlives_g,IPD_Control%nfracseed,&
IPD_Control%nseed_g,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1),IPD_Control%nsmooth,IPD_Control%ca_amplitude)
endif
endif

!--- if coupled, assign coupled fields
Expand Down Expand Up @@ -656,12 +660,21 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
call mpp_error(FATAL, 'Logic errror: cellular_automata not compatible with non-uniform blocksizes')
end if
! *DH
call cellular_automata(IPD_Control%kdt, IPD_Data(:)%Statein, IPD_Data(:)%Coupling, IPD_Data(:)%Intdiag, &
Atm_block%nblks, IPD_Control%levs, IPD_Control%nca, IPD_Control%ncells, &
IPD_Control%nlives, IPD_Control%nfracseed, IPD_Control%nseed, &
IPD_Control%nthresh, IPD_Control%ca_global, IPD_Control%ca_sgs, &
IPD_Control%iseed_ca, IPD_Control%ca_smooth, IPD_Control%nspinup, &
Atm_block%blksz(1))
if(IPD_Control%do_ca)then
if(IPD_Control%ca_sgs == .true.)then
call cellular_automata_sgs(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
IPD_Control%nca,IPD_Control%ncells,IPD_Control%nlives,IPD_Control%nfracseed,&
IPD_Control%nseed,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1))
endif
if(IPD_Control%ca_global == .true.)then
call cellular_automata_global(IPD_Control%kdt,IPD_Data(:)%Statein,IPD_Data(:)%Coupling,IPD_Data(:)%Intdiag,Atm_block%nblks,IPD_Control%levs, &
IPD_Control%nca_g,IPD_Control%ncells_g,IPD_Control%nlives_g,IPD_Control%nfracseed,&
IPD_Control%nseed_g,IPD_Control%nthresh,IPD_Control%ca_global,IPD_Control%ca_sgs,IPD_Control%iseed_ca,&
IPD_Control%ca_smooth,IPD_Control%nspinup,Atm_block%blksz(1),IPD_Control%nsmooth,IPD_Control%ca_amplitude)
endif

endif
endif

Atm(mytile)%flagstruct%do_skeb = IPD_Control%do_skeb
Expand Down
2 changes: 1 addition & 1 deletion ccpp/physics
Submodule physics updated 11 files
+9 −37 physics/GFS_DCNV_generic.F90
+0 −77 physics/GFS_DCNV_generic.meta
+3 −3 physics/GFS_MP_generic.F90
+8 −0 physics/GFS_MP_generic.meta
+1 −6 physics/GFS_debug.F90
+78 −34 physics/GFS_stochastics.F90
+51 −0 physics/GFS_stochastics.meta
+3 −3 physics/GFS_surface_generic.F90
+8 −0 physics/GFS_surface_generic.meta
+82 −14 physics/samfdeepcnv.f
+42 −0 physics/samfdeepcnv.meta
4 changes: 2 additions & 2 deletions gfsphysics/GFS_layer/GFS_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2028,13 +2028,13 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop

idx = idx + 1
ExtDiag(idx)%axes = 2
ExtDiag(idx)%name = 'ca_out'
ExtDiag(idx)%name = 'ca1'
ExtDiag(idx)%desc = 'Cellular Automata'
ExtDiag(idx)%unit = '%'
ExtDiag(idx)%mod_name = 'gfs_phys'
allocate (ExtDiag(idx)%data(nblks))
do nb = 1,nblks
ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ca_out(:)
ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ca1(:)
enddo

idx = idx + 1
Expand Down
128 changes: 89 additions & 39 deletions gfsphysics/GFS_layer/GFS_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, &
integer :: ntrac
integer :: ix
#ifndef CCPP
integer :: blocksize
integer :: blocksize,k
real(kind=kind_phys), allocatable :: si(:)
real(kind=kind_phys), parameter :: p_ref = 101325.0d0
#endif
Expand Down Expand Up @@ -194,10 +194,9 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, &
Init_parm%iau_offset, Init_parm%bdat, &
Init_parm%cdat, Init_parm%tracer_names, &
Init_parm%input_nml_file, Init_parm%tile_num, &
Init_parm%blksz &
Init_parm%blksz,Init_parm%ak, Init_parm%bk &
#ifdef CCPP
,Init_parm%ak, Init_parm%bk, &
Init_parm%restart, Init_parm%hydrostatic, &
,Init_parm%restart, Init_parm%hydrostatic, &
communicator, ntasks, nthrds &
#endif
)
Expand Down Expand Up @@ -469,6 +468,29 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, &
!--- this note is placed here to alert users to study
!--- the FV3GFS_io.F90 module

#ifndef CCPP
if(Model%do_ca .and. Model%ca_global)then

do nb = 1,nblks
do k=1,Model%levs
if (Model%si(k) .lt. 0.1 .and. Model%si(k) .gt. 0.025) then
Coupling(nb)%vfact_ca(k) = (Model%si(k)-0.025)/(0.1-0.025)
else if (Model%si(k) .lt. 0.025) then
Coupling(nb)%vfact_ca(k) = 0.0
else
Coupling(nb)%vfact_ca(k) = 1.0
endif
enddo
enddo

do nb = 1,nblks
Coupling(nb)%vfact_ca(2)=Coupling(nb)%vfact_ca(3)*0.5
Coupling(nb)%vfact_ca(1)=0.0
enddo

endif
#endif

end subroutine GFS_initialize


Expand Down Expand Up @@ -710,8 +732,7 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
!--- local variables
integer :: k, i
real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew,sppt_vwt
!real(kind=kind_phys),dimension(size(Statein%tgrs,1),size(Statein%tgrs,2)) :: tconvtend, &
! qconvtend,uconvtend,vconvtend
real(kind=kind_phys),dimension(size(Statein%tgrs,1),size(Statein%tgrs,2)) :: ca1

if (Model%do_sppt) then
do k = 1,size(Statein%tgrs,2)
Expand Down Expand Up @@ -739,27 +760,11 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
Diag%sppt_wts(i,Model%levs-k+1)=Coupling%sppt_wts(i,k)


! if(Model%isppt_deep)then

! tconvtend(i,k)=Coupling%tconvtend(i,k)
! qconvtend(i,k)=Coupling%qconvtend(i,k)
! uconvtend(i,k)=Coupling%uconvtend(i,k)
! vconvtend(i,k)=Coupling%vconvtend(i,k)


! upert = (Stateout%gu0(i,k) - Statein%ugrs(i,k) - uconvtend(i,k)) + uconvtend(i,k) * Coupling%sppt_wts(i,k)
! vpert = (Stateout%gv0(i,k) - Statein%vgrs(i,k) - vconvtend(i,k)) + vconvtend(i,k) * Coupling%sppt_wts(i,k)
! tpert = (Stateout%gt0(i,k) - Statein%tgrs(i,k) - Tbd%dtdtr(i,k) - tconvtend(i,k)) + tconvtend(i,k) * Coupling%sppt_wts(i,k)
! qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1) - qconvtend(i,k)) + qconvtend(i,k) * Coupling%sppt_wts(i,k)

!else

upert = (Stateout%gu0(i,k) - Statein%ugrs(i,k)) * Coupling%sppt_wts(i,k)
vpert = (Stateout%gv0(i,k) - Statein%vgrs(i,k)) * Coupling%sppt_wts(i,k)
tpert = (Stateout%gt0(i,k) - Statein%tgrs(i,k) - Tbd%dtdtr(i,k)) * Coupling%sppt_wts(i,k)
qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1)) * Coupling%sppt_wts(i,k)

!endif

Stateout%gu0(i,k) = Statein%ugrs(i,k)+upert
Stateout%gv0(i,k) = Statein%vgrs(i,k)+vpert
Expand All @@ -773,22 +778,6 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
enddo
enddo

!if(Model%isppt_deep == .true.)then
! Sfcprop%tprcp(:) = Sfcprop%tprcp(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
! Diag%totprcp(:) = Diag%totprcp(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
! Diag%cnvprcp(:) = Diag%cnvprcp(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
!! ! bucket precipitation adjustment due to sppt
! Diag%totprcpb(:) = Diag%totprcpb(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)
! Diag%cnvprcpb(:) = Diag%cnvprcpb(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rainc(:)


! if (Model%cplflx) then !Need to make proper adjustments for deep convection only perturbations
! Coupling%rain_cpl(:) = Coupling%rain_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%drain_cpl(:)
! Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%dsnow_cpl(:)
! endif

!else

! instantaneous precip rate going into land model at the next time step
Sfcprop%tprcp(:) = Coupling%sppt_wts(:,15)*Sfcprop%tprcp(:)
Diag%totprcp(:) = Diag%totprcp(:) + (Coupling%sppt_wts(:,15) - 1 )*Diag%rain(:)
Expand All @@ -804,10 +793,71 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, &
Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (Coupling%sppt_wts(:,15) - 1.0)*Tbd%dsnow_cpl(:)
endif

!endif
endif


if (Model%do_ca .and. Model%ca_global) then
do k = 1,size(Statein%tgrs,2)
do i = 1,size(Statein%tgrs,1)
sppt_vwt=1.0
if (Diag%zmtnblck(i).EQ.0.0) then
sppt_vwt=1.0
else
if (k.GT.Diag%zmtnblck(i)+2) then
sppt_vwt=1.0
endif
if (k.LE.Diag%zmtnblck(i)) then
sppt_vwt=0.0
endif
if (k.EQ.Diag%zmtnblck(i)+1) then
sppt_vwt=0.333333
endif
if (k.EQ.Diag%zmtnblck(i)+2) then
sppt_vwt=0.666667
endif
endif

ca1(i,k)=((Coupling%ca1(i)-1.)*sppt_vwt*Coupling%vfact_ca(k))+1.0

upert = (Stateout%gu0(i,k) - Statein%ugrs(i,k)) * ca1(i,k)
vpert = (Stateout%gv0(i,k) - Statein%vgrs(i,k)) * ca1(i,k)
tpert = (Stateout%gt0(i,k) - Statein%tgrs(i,k) - Tbd%dtdtr(i,k)) * ca1(i,k)
qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1)) * ca1(i,k)

Stateout%gu0(i,k) = Statein%ugrs(i,k)+upert
Stateout%gv0(i,k) = Statein%vgrs(i,k)+vpert

!negative humidity check
qnew = Statein%qgrs(i,k,1)+qpert
if (qnew >= 1.0e-10) then
Stateout%gq0(i,k,1) = qnew
Stateout%gt0(i,k) = Statein%tgrs(i,k) + tpert + Tbd%dtdtr(i,k)
endif

enddo
enddo



! instantaneous precip rate going into land model at the next time step
Sfcprop%tprcp(:) = ca1(:,15)*Sfcprop%tprcp(:)
Diag%totprcp(:) = Diag%totprcp(:) + (ca1(:,15) - 1 )*Diag%rain(:)
! acccumulated total and convective preciptiation
Diag%cnvprcp(:) = Diag%cnvprcp(:) + (ca1(:,15) - 1 )*Diag%rainc(:)
! bucket precipitation adjustment due to sppt
Diag%totprcpb(:) = Diag%totprcpb(:) + (ca1(:,15) - 1 )*Diag%rain(:)
Diag%cnvprcpb(:) = Diag%cnvprcpb(:) + (ca1(:,15) - 1 )*Diag%rainc(:)

if (Model%cplflx) then
Coupling%rain_cpl(:) = Coupling%rain_cpl(:) + (ca1(:,15) - 1.0)*Tbd%drain_cpl(:)
Coupling%snow_cpl(:) = Coupling%snow_cpl(:) + (ca1(:,15) - 1.0)*Tbd%dsnow_cpl(:)
endif

endif




if (Model%do_shum) then
Stateout%gq0(:,:,1) = Stateout%gq0(:,:,1)*(1.0 + Coupling%shum_wts(:,:))
endif
Expand Down
Loading

0 comments on commit e8fd6ee

Please sign in to comment.