Skip to content
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.

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

upgraded OPAC aerosol with more advanced MERRA2 aerosol #594

Merged
merged 7 commits into from
Mar 26, 2021
Merged
29 changes: 21 additions & 8 deletions physics/GFS_phys_time_vary.fv3.F90
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ module GFS_phys_time_vary
use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol

use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm
use aerinterp, only : read_aerdata, setindxaer, aerinterpol
use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf

use iccn_def, only : ciplin, ccnin, ci_pres
use iccninterp, only : read_cidata, setindxci, ciinterpol
Expand Down Expand Up @@ -166,7 +166,7 @@ subroutine GFS_phys_time_vary_init (
integer, intent(out) :: errflg

! Local variables
integer :: i, j, ix, vegtyp
integer :: i, j, ix, vegtyp, iamin, iamax, jamin, jamax
real(kind_phys) :: rsnow

!--- Noah MP
Expand All @@ -182,12 +182,17 @@ subroutine GFS_phys_time_vary_init (
errflg = 0

if (is_initialized) return
iamin=999
climbfuji marked this conversation as resolved.
Show resolved Hide resolved
iamax=-999
jamin=999
jamax=-999

!$OMP parallel num_threads(nthrds) default(none) &
!$OMP shared (me,master,ntoz,h2o_phys,im,nx,ny,idate) &
!$OMP shared (xlat_d,xlon_d,imap,jmap,errmsg,errflg) &
!$OMP shared (levozp,oz_coeff,oz_pres,ozpl) &
!$OMP shared (levh2o,h2o_coeff,h2o_pres,h2opl) &
!$OMP shared (iamin, iamax, jamin, jamax) &
!$OMP shared (iaerclm,ntrcaer,aer_nm,iflip,iccn) &
!$OMP shared (jindx1_o3,jindx2_o3,ddy_o3,jindx1_h,jindx2_h,ddy_h) &
!$OMP shared (jindx1_aer,jindx2_aer,ddy_aer,iindx1_aer,iindx2_aer,ddx_aer) &
Expand Down Expand Up @@ -300,16 +305,20 @@ subroutine GFS_phys_time_vary_init (
call setindxh2o (im, xlat_d, jindx1_h, jindx2_h, ddy_h)
endif

!$OMP section
!> - Call setindxaer() to initialize aerosols data
!$OMP section
if (iaerclm) then
call setindxaer (im, xlat_d, jindx1_aer, &
jindx2_aer, ddy_aer, xlon_d, &
iindx1_aer, iindx2_aer, ddx_aer, &
me, master)
iamin=min(minval(iindx1_aer), iamin)
iamax=max(maxval(iindx2_aer), iamax)
jamin=min(minval(jindx1_aer), jamin)
jamax=max(maxval(jindx2_aer), jamax)
endif

!$OMP section

!> - Call setindxci() to initialize IN and CCN data
if (iccn == 1) then
call setindxci (im, xlat_d, jindx1_ci, &
Expand Down Expand Up @@ -367,6 +376,10 @@ subroutine GFS_phys_time_vary_init (
!$OMP end sections

!$OMP end parallel
if (iaerclm) then
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is ok to leave it here for the moment, given that I will be working on speeding up the code in physics/GFS_phys_time_vary.fv3.F90 after your commit went in.

call read_aerdataf (iamin, iamax, jamin, jamax, me,master,iflip, &
idate,errmsg,errflg)
endif

if (lsm == lsm_noahmp) then
if (all(tvxy < zero)) then
Expand Down Expand Up @@ -804,10 +817,10 @@ subroutine GFS_phys_time_vary_timestep_init (

!> - Call ciinterpol() to make IN and CCN data interpolation
if (iccn == 1) then
call ciinterpol (me, im, idate, fhour, &
jindx1_ci, jindx2_ci, &
ddy_ci, iindx1_ci, &
iindx2_ci, ddx_ci, &
call ciinterpol (me, im, idate, fhour, &
jindx1_ci, jindx2_ci, &
ddy_ci, iindx1_ci, &
iindx2_ci, ddx_ci, &
levs, prsl, in_nm, ccn_nm)
endif

Expand Down
4 changes: 2 additions & 2 deletions physics/aerclm_def.F
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@ module aerclm_def
use machine , only : kind_phys
implicit none

integer, parameter :: levsaer=50, ntrcaerm=15, timeaer=12
integer :: latsaer, lonsaer, ntrcaer
integer, parameter :: levsaer=72, ntrcaerm=15, timeaer=12
integer :: latsaer, lonsaer, ntrcaer, levsw

character*10 :: specname(ntrcaerm)
real (kind=kind_phys):: aer_time(13)
Expand Down
99 changes: 48 additions & 51 deletions physics/aerinterp.F90
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ module aerinterp

private

public :: read_aerdata, setindxaer, aerinterpol
public :: read_aerdata, setindxaer, aerinterpol, read_aerdataf

contains

Expand All @@ -32,11 +32,6 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
logical :: file_exist

integer, allocatable :: invardims(:)
real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff
real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx
real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp
real(kind=kind_io8),allocatable,dimension(:) :: aer_lati
real(kind=kind_io8),allocatable,dimension(:) :: aer_loni
!
!! ===================================================================
if (me == master) then
Expand Down Expand Up @@ -72,50 +67,62 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
! specify latsaer, lonsaer, hmx
lonsaer = dim1
latsaer = dim2
hmx = int(dim1/2) ! to swap long from W-E to E-W
levsw = dim3

if(me==master) then
print *, 'MERRA2 dim: ',dim1, dim2, dim3
endif

! allocate arrays
if (.not. allocated(aer_loni)) then
allocate (aer_loni(lonsaer))
allocate (aer_lati(latsaer))
endif

if (.not. allocated(aer_lat)) then
allocate(aer_lat(latsaer))
allocate(aer_lon(lonsaer))
allocate(aerin(lonsaer,latsaer,levsaer,ntrcaerm,timeaer))
allocate(aer_pres(lonsaer,latsaer,levsaer,timeaer))
endif

! construct lat/lon array
call nf_inq_varid(ncid, 'lat', varid)
call nf_get_var(ncid, varid, aer_lati)
call nf_get_var(ncid, varid, aer_lat)
call nf_inq_varid(ncid, 'lon', varid)
call nf_get_var(ncid, varid, aer_loni)

do i = 1, hmx ! flip from (-180,180) to (0,360)
if(aer_loni(i)<0.) aer_loni(i)=aer_loni(i)+360.
aer_lon(i+hmx) = aer_loni(i)
aer_lon(i) = aer_loni(i+hmx)
enddo
call nf_get_var(ncid, varid, aer_lon)
call nf_close(ncid)
END SUBROUTINE read_aerdata
!
!**********************************************************************
SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, &
me, master, iflip, idate, errmsg, errflg)
use machine, only: kind_phys, kind_io4, kind_io8
use aerclm_def
use netcdf

do i = 1, latsaer
aer_lat(i) = aer_lati(i)
enddo
!--- in/out
integer, intent(in) :: me, master, iflip, idate(4)
integer, intent(in) :: iamin, iamax, jamin, jamax
character(len=*), intent(inout) :: errmsg
integer, intent(inout) :: errflg

call nf_close(ncid)
!--- locals
integer :: ncid, varid
integer :: i, j, k, n, ii, imon, klev
character :: fname*50, mn*2, vname*10
logical :: file_exist
integer, allocatable :: invardims(:)
real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff
real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx
real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp
!
if (.not. allocated(aerin)) then
allocate(aerin(iamin:iamax,jamin:jamax,levsaer,ntrcaerm,timeaer))
allocate(aer_pres(iamin:iamax,jamin:jamax,levsaer,timeaer))
endif

! allocate local working arrays
if (.not. allocated(buff)) then
allocate (buff(lonsaer, latsaer, dim3))
allocate (pres_tmp(lonsaer,dim3))
allocate (buff(lonsaer, latsaer, levsw))
allocate (pres_tmp(lonsaer,levsw))
endif
if (.not. allocated(buffx)) then
allocate (buffx(lonsaer, latsaer, dim3,1))
allocate (buffx(lonsaer, latsaer, levsw,1))
endif

!! ===================================================================
Expand All @@ -137,11 +144,11 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
call nf_inq_varid(ncid, "DELP", varid)
call nf_get_var(ncid, varid, buff)

do j = 1, latsaer
do i = 1, lonsaer
do j = jamin, jamax
do i = iamin, iamax
! constract pres_tmp (top-down), note input is top-down
pres_tmp(i,1) = 0.
do k=2, dim3
do k=2, levsw
pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k)
enddo !k-loop
enddo !i-loop (lon)
Expand All @@ -151,11 +158,10 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
if ( iflip == 0 ) then ! data from toa to sfc
klev = k
else ! data from sfc to top
klev = ( dim3 - k ) + 1
klev = ( levsw - k ) + 1
endif
do i = 1, hmx
aer_pres(i+hmx,j,k,imon)= 1.d0*pres_tmp(i,klev)
aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i+hmx,klev)
do i = iamin, iamax
aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i,klev)
enddo !i-loop (lon)
enddo !k-loop (lev)
enddo !j-loop (lat)
Expand All @@ -168,22 +174,18 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
call nf_inq_varid(ncid, vname, varid)
call nf_get_var(ncid, varid, buffx)

do j = 1, latsaer
do j = jamin, jamax
do k = 1, levsaer
! input is from toa to sfc
if ( iflip == 0 ) then ! data from toa to sfc
klev = k
else ! data from sfc to top
klev = ( dim3 - k ) + 1
klev = ( levsw - k ) + 1
endif
do i = 1, hmx
aerin(i+hmx,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1)
if(aerin(i+hmx,j,k,ii,imon)<0.or.aerin(i+hmx,j,k,ii,imon)>1.) then
aerin(i+hmx,j,k,ii,imon) = 0.
end if
aerin(i,j,k,ii,imon) = 1.d0*buffx(i+hmx,j,klev,1)
do i = iamin, iamax
aerin(i,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1)
if(aerin(i,j,k,ii,imon)<0.or.aerin(i,j,k,ii,imon)>1.) then
aerin(i,j,k,ii,imon) = 0.
aerin(i,j,k,ii,imon) = 1.e-15
end if
enddo !i-loop (lon)
enddo !k-loop (lev)
Expand All @@ -195,13 +197,9 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg)
call nf_close(ncid)
enddo !imon-loop
!---
deallocate (aer_loni, aer_lati)
deallocate (buff, pres_tmp)
deallocate (buffx)

END SUBROUTINE read_aerdata
!
!**********************************************************************
END SUBROUTINE read_aerdataf
!
SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, &
iindx1,iindx2,ddx,me,master)
Expand Down Expand Up @@ -341,7 +339,6 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, &
+TEMI*DDY(j)*aer_pres(I1,J2,L,n1)+DDX(j)*TEMJ*aer_pres(I2,J1,L,n1))&
+tx2*(TEMI*TEMJ*aer_pres(I1,J1,L,n2)+DDX(j)*DDY(J)*aer_pres(I2,J2,L,n2) &
+TEMI*DDY(j)*aer_pres(I1,J2,L,n2)+DDX(j)*TEMJ*aer_pres(I2,J1,L,n2))

ENDDO
ENDDO

Expand Down Expand Up @@ -369,7 +366,7 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, &
tx1 = temi/(aerpres(j,i1) - aerpres(j,i2))
tx2 = temj/(aerpres(j,i1) - aerpres(j,i2))
DO ii = 1, ntrcaer
aerout(j,L,ii)= aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2
aerout(j,L,ii)= aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2
ENDDO
endif
ENDDO !L-loop
Expand Down
4 changes: 3 additions & 1 deletion physics/radiation_aerosols.f
Original file line number Diff line number Diff line change
Expand Up @@ -561,7 +561,7 @@ subroutine aer_init &

laswflg= (mod(iaerflg,10) > 0) ! control flag for sw tropospheric aerosol
lalwflg= (mod(iaerflg/10,10) > 0) ! control flag for lw tropospheric aerosol
lavoflg= (iaerflg >= 100) ! control flag for stratospheric volcanic aeros
lavoflg= (mod(iaerflg/100,10) >0) ! control flag for stratospheric volcanic aeros

!> -# Call wrt_aerlog() to write aerosol parameter configuration to output logs.

Expand Down Expand Up @@ -4446,6 +4446,8 @@ subroutine aeropt
asy1 = f_zero
sca1 = f_zero
ssa1 = f_zero
asy = f_zero
ssa = f_zero
do m = 1, kcm1
cm = max(aerms(k,m),0.0) * dz1(k)
ext1 = ext1 + cm*extrhi_grt(m,ib)
Expand Down
8 changes: 4 additions & 4 deletions physics/sfc_sice.f
Original file line number Diff line number Diff line change
Expand Up @@ -287,11 +287,11 @@ subroutine sfc_sice_run &
q0 = min(qs1, q0)

if (fice(i) < cimin) then
print *,'warning: ice fraction is low:', fice(i)
! print *,'warning: ice fraction is low:', fice(i)
fice(i) = cimin
tice(i) = tgice
tskin(i)= tgice
print *,'fix ice fraction: reset it to:', fice(i)
! print *,'fix ice fraction: reset it to:', fice(i)
endif
ffw(i) = one - fice(i)

Expand Down Expand Up @@ -362,9 +362,9 @@ subroutine sfc_sice_run &
snowd(i) = min( snowd(i), hsmax )

if (snowd(i) > (2.0_kind_phys*hice(i))) then
print *, 'warning: too much snow :',snowd(i)
! print *, 'warning: too much snow :',snowd(i)
snowd(i) = hice(i) + hice(i)
print *,'fix: decrease snow depth to:',snowd(i)
! print *,'fix: decrease snow depth to:',snowd(i)
endif
endif
enddo
Expand Down