Skip to content

Commit

Permalink
Merge pull request #594 from AnningCheng-NOAA/merra2
Browse files Browse the repository at this point in the history 
upgraded OPAC aerosol with more advanced MERRA2 aerosol
  • Loading branch information
@climbfuji
climbfuji authored Mar 26, 2021
2 parents 0f9b52d + b6ef6c3 commit affd4a6
Show file tree
Hide file tree
Showing 5 changed files with 78 additions and 66 deletions.
29 changes: 21 additions & 8 deletions physics/GFS_phys_time_vary.fv3.F90
View file
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
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
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit affd4a6

Please sign in to comment.