From 9e784fd219a60457130995a8e94fbe582b10153f Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 9 Nov 2021 09:47:34 -0500 Subject: [PATCH 1/6] decreases memory usage by 6 times for MERRA2 --- physics/GFS_phys_time_vary.fv3.F90 | 18 +-- physics/GFS_phys_time_vary.fv3.meta | 15 ++ physics/aerclm_def.F | 4 +- physics/aerinterp.F90 | 241 ++++++++++++++++++++++++---- 4 files changed, 241 insertions(+), 37 deletions(-) diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index a8ecc1a5e..d6155e6b1 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -21,7 +21,7 @@ module GFS_phys_time_vary use h2o_def, only : levh2o, h2o_coeff, h2o_lat, h2o_pres, h2o_time, h2oplin use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol - use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm + use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm, iamin, iamax, jamin, jamax use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf use iccn_def, only : ciplin, ccnin, ci_pres @@ -68,7 +68,7 @@ module GFS_phys_time_vary !! @{ subroutine GFS_phys_time_vary_init ( & me, master, ntoz, h2o_phys, iaerclm, iccn, iflip, im, nx, ny, idate, xlat_d, xlon_d, & - jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & + jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl,fhour, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap, & do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau, & @@ -88,6 +88,7 @@ subroutine GFS_phys_time_vary_init ( integer, intent(in) :: me, master, ntoz, iccn, iflip, im, nx, ny logical, intent(in) :: h2o_phys, iaerclm, flag_restart integer, intent(in) :: idate(:) + real(kind_phys), intent(in) :: fhour real(kind_phys), intent(in) :: xlat_d(:), xlon_d(:) integer, intent(inout) :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:) @@ -173,7 +174,7 @@ subroutine GFS_phys_time_vary_init ( integer, intent(out) :: errflg ! Local variables - integer :: i, j, ix, vegtyp, iamin, iamax, jamin, jamax + integer :: i, j, ix, vegtyp real(kind_phys) :: rsnow !--- Noah MP @@ -387,8 +388,7 @@ subroutine GFS_phys_time_vary_init ( if (errflg/=0) return if (iaerclm) then - call read_aerdataf (iamin, iamax, jamin, jamax, me, master, iflip, & - idate, errmsg, errflg) + call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg) if (errflg/=0) return end if @@ -715,7 +715,7 @@ end subroutine GFS_phys_time_vary_init subroutine GFS_phys_time_vary_timestep_init ( & me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, nsswr, fhswr, lsswr, fhour, & imfdeepcnv, cal_pre, random_clds, nscyc, ntoz, h2o_phys, iaerclm, iccn, clstp, & - jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & + jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, iflip, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, & imap, jmap, prsl, seed0, rann, nthrds, nx, ny, nsst, tile_num, nlunit, lsoil, lsoil_lsm,& @@ -730,7 +730,7 @@ subroutine GFS_phys_time_vary_timestep_init ( ! Interface variables integer, intent(in) :: me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, & - nsswr, imfdeepcnv, iccn, nscyc, ntoz + nsswr, imfdeepcnv, iccn, nscyc, ntoz, iflip integer, intent(in) :: idate(:) real(kind_phys), intent(in) :: fhswr, fhour logical, intent(in) :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm @@ -797,7 +797,7 @@ subroutine GFS_phys_time_vary_timestep_init ( !$OMP shared(ozpl,ddy_o3,h2o_phys,jindx1_h,jindx2_h,h2opl,ddy_h,iaerclm,master) & !$OMP shared(levs,prsl,iccn,jindx1_ci,jindx2_ci,ddy_ci,iindx1_ci,iindx2_ci) & !$OMP shared(ddx_ci,in_nm,ccn_nm,do_ugwp_v1,jindx1_tau,jindx2_tau,ddy_j1tau) & -!$OMP shared(ddy_j2tau,tau_amf) & +!$OMP shared(ddy_j2tau,tau_amf,iflip) & !$OMP private(iseed,iskip,i,j,k) !$OMP sections @@ -889,7 +889,7 @@ subroutine GFS_phys_time_vary_timestep_init ( ! aerinterpol is using threading inside, don't ! move into OpenMP parallel section above call aerinterpol (me, master, nthrds, im, idate, & - fhour, jindx1_aer, jindx2_aer, & + fhour, iflip, jindx1_aer, jindx2_aer, & ddy_aer, iindx1_aer, & iindx2_aer, ddx_aer, & levs, prsl, aer_nm) diff --git a/physics/GFS_phys_time_vary.fv3.meta b/physics/GFS_phys_time_vary.fv3.meta index 6c7f086dd..90f061dd1 100644 --- a/physics/GFS_phys_time_vary.fv3.meta +++ b/physics/GFS_phys_time_vary.fv3.meta @@ -162,6 +162,14 @@ type = real kind = kind_phys intent = in +[fhour] + standard_name = forecast_time + long_name = current forecast time + units = h + dimensions = () + type = real + kind = kind_phys + intent = in [jindx1_aer] standard_name = lower_latitude_index_of_aerosol_forcing_for_interpolation long_name = interpolation low index for prescribed aerosols in the y direction @@ -1151,6 +1159,13 @@ type = real kind = kind_phys intent = inout +[iflip] + standard_name = control_for_vertical_index_direction + long_name = iflip - is not the same as flipv + units = flag + dimensions = () + type = integer + intent = in [jindx1_aer] standard_name = lower_latitude_index_of_aerosol_forcing_for_interpolation long_name = interpolation low index for prescribed aerosols in the y direction diff --git a/physics/aerclm_def.F b/physics/aerclm_def.F index e66825278..b6760f30c 100644 --- a/physics/aerclm_def.F +++ b/physics/aerclm_def.F @@ -2,8 +2,10 @@ module aerclm_def use machine , only : kind_phys, kind_io4 implicit none - integer, parameter :: levsaer=72, ntrcaerm=15, timeaer=12 + integer, parameter :: levsaer=72, ntrcaerm=15, timeaer=2 integer :: latsaer, lonsaer, ntrcaer, levsw + integer :: n1sv, n2sv + integer :: iamin, iamax, jamin, jamax character*10 :: specname(ntrcaerm) real (kind=kind_phys):: aer_time(13) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 4b3232ab1..a99fae7b0 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -11,7 +11,7 @@ module aerinterp private - public :: read_aerdata, setindxaer, aerinterpol, read_aerdataf + public :: read_aerdata, setindxaer, aerinterpol,read_aerdataf contains @@ -45,7 +45,7 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) !! =================================================================== !! check if all necessary files exist !! =================================================================== - do imon = 1, timeaer + do imon = 1, 12 write(mn,'(i2.2)') imon fname=trim("aeroclim.m"//mn//".nc") inquire (file = fname, exist = file_exist) @@ -97,23 +97,27 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) END SUBROUTINE read_aerdata ! !********************************************************************** - SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & - me, master, iflip, idate, errmsg, errflg) + SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) use machine, only: kind_phys, kind_io4, kind_io8 use aerclm_def use netcdf !--- 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 - + real(kind=kind_phys), intent(in) :: fhour !--- locals integer :: ncid, varid - integer :: i, j, k, n, ii, imon, klev + integer :: i, j, k, n, ii, imon, klev, n1, n2 character :: fname*50, mn*2, vname*10 logical :: file_exist + integer IDAT(8),JDAT(8) + real(kind=kind_phys) RINC(5), rjday + integer jdow, jdoy, jday + real(4) rinc4(5) + integer w3kindreal,w3kindint + integer, allocatable :: invardims(:) real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx @@ -128,12 +132,104 @@ SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & allocate (buff(lonsaer, latsaer, levsw)) allocate (pres_tmp(lonsaer, levsw)) allocate (buffx(lonsaer, latsaer, levsw, 1)) - +!! found interpolation months + IDAT = 0 + IDAT(1) = IDATE(4) + IDAT(2) = IDATE(2) + IDAT(3) = IDATE(3) + IDAT(5) = IDATE(1) + RINC = 0. + RINC(2) = FHOUR + call w3kind(w3kindreal,w3kindint) + if(w3kindreal == 4) then + rinc4 = rinc + CALL W3MOVDAT(RINC4,IDAT,JDAT) + else + CALL W3MOVDAT(RINC,IDAT,JDAT) + endif +! + jdow = 0 + jdoy = 0 + jday = 0 + call w3doxdat(jdat,jdow,jdoy,jday) + rjday = jdoy + jdat(5) / 24. + IF (RJDAY < aer_time(1)) RJDAY = RJDAY+365. +! + n2 = 13 + do j=2, 12 + if (rjday < aer_time(j)) then + n2 = j + exit + endif + enddo + n1 = n2 - 1 + if (n2 > 12) n2 = n2 -12 + write(*,*)"AAA0",n1, n2, iamin, iamax, jamin, jamax !! =================================================================== !! loop thru m01 - m12 for aer/pres array !! =================================================================== - do imon = 1, timeaer - write(mn,'(i2.2)') imon + write(mn,'(i2.2)') n1 + fname=trim("aeroclim.m"//mn//".nc") + call nf_open(fname , nf90_NOWRITE, ncid) + +! ====> construct 3-d pressure array (Pa) + call nf_inq_varid(ncid, "DELP", varid) + call nf_get_var(ncid, varid, buff) + + 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, levsw + pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) + enddo !k-loop + enddo !i-loop (lon) + +! extract pres_tmp to fill aer_pres (in Pa) + do k = 1, levsaer + if ( iflip == 0 ) then ! data from toa to sfc + klev = k + else ! data from sfc to top + klev = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aer_pres(i,j,k,1) = 1.d0*pres_tmp(i,klev) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + +! ====> construct 4-d aerosol array (kg/kg) +! merra2 data is top down +! for GFS, iflip 0: toa to sfc; 1: sfc to toa + DO ii = 1, ntrcaerm + vname=trim(specname(ii)) + call nf_inq_varid(ncid, vname, varid) + write(*,*)"AAA2",vname + call nf_get_var(ncid, varid, buffx) + + 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 = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aerin(i,j,k,ii,1) = 1.d0*buffx(i,j,klev,1) + if(aerin(i,j,k,ii,1) < 0 .or. aerin(i,j,k,ii,1) > 1.) then + aerin(i,j,k,ii,1) = 1.e-15 + endif + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + + ENDDO ! ii-loop (ntracaerm) + +! close the file + call nf_close(ncid) +!! =================================================================== + write(mn,'(i2.2)') n2 fname=trim("aeroclim.m"//mn//".nc") call nf_open(fname , nf90_NOWRITE, ncid) @@ -158,7 +254,7 @@ SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & klev = ( levsw - k ) + 1 endif do i = iamin, iamax - aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i,klev) + aer_pres(i,j,k,2) = 1.d0*pres_tmp(i,klev) enddo !i-loop (lon) enddo !k-loop (lev) enddo !j-loop (lat) @@ -180,9 +276,9 @@ SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & klev = ( levsw - k ) + 1 endif 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) = 1.e-15 + aerin(i,j,k,ii,2) = 1.d0*buffx(i,j,klev,1) + if(aerin(i,j,k,ii,2) < 0 .or. aerin(i,j,k,ii,2) > 1.) then + aerin(i,j,k,ii,2) = 1.e-15 endif enddo !i-loop (lon) enddo !k-loop (lev) @@ -192,7 +288,8 @@ SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & ! close the file call nf_close(ncid) - enddo !imon-loop + n1sv=n1 + n2sv=n2 !--- deallocate (buff, pres_tmp) deallocate (buffx) @@ -256,15 +353,20 @@ END SUBROUTINE setindxaer !********************************************************************** !********************************************************************** ! - SUBROUTINE aerinterpol(me,master,nthrds,npts,IDATE,FHOUR,jindx1,jindx2, & + SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, & ddy,iindx1,iindx2,ddx,lev,prsl,aerout) ! - USE MACHINE, ONLY : kind_phys + use machine, only: kind_phys, kind_io4, kind_io8 use aerclm_def + use netcdf + implicit none - integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii + integer, intent(in) :: iflip + integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii, klev real(kind=kind_phys) fhour,temj, tx1, tx2,temi, tem real(kind=kind_phys), dimension(npts) :: temij,temiy,temjx,ddxy + character :: fname*50, mn*2, vname*10 + ! integer JINDX1(npts), JINDX2(npts), iINDX1(npts), iINDX2(npts) @@ -279,6 +381,11 @@ SUBROUTINE aerinterpol(me,master,nthrds,npts,IDATE,FHOUR,jindx1,jindx2, & integer jdow, jdoy, jday real(4) rinc4(5) integer w3kindreal,w3kindint + integer ncid, varid + real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff + real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx + real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp + ! IDAT = 0 IDAT(1) = IDATE(4) @@ -310,6 +417,86 @@ SUBROUTINE aerinterpol(me,master,nthrds,npts,IDATE,FHOUR,jindx1,jindx2, & endif enddo n1 = n2 - 1 +! need to read a new month + if (n1>n1sv) then + DO ii = 1, ntrcaerm + do j = jamin, jamax + do k = 1, levsaer + do i = iamin, iamax + aerin(i,j,k,ii,1) = aerin(i,j,k,ii,2) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + ENDDO ! ii-loop (ntracaerm) +!! =================================================================== + allocate (buff(lonsaer, latsaer, levsw)) + allocate (pres_tmp(lonsaer, levsw)) + allocate (buffx(lonsaer, latsaer, levsw, 1)) + + write(mn,'(i2.2)') n2 + fname=trim("aeroclim.m"//mn//".nc") + call nf_open(fname , nf90_NOWRITE, ncid) + +! ====> construct 3-d pressure array (Pa) + call nf_inq_varid(ncid, "DELP", varid) + call nf_get_var(ncid, varid, buff) + + 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, levsw + pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) + enddo !k-loop + enddo !i-loop (lon) + +! extract pres_tmp to fill aer_pres (in Pa) + do k = 1, levsaer + if ( iflip == 0 ) then ! data from toa to sfc + klev = k + else ! data from sfc to top + klev = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aer_pres(i,j,k,2) = 1.d0*pres_tmp(i,klev) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + +! ====> construct 4-d aerosol array (kg/kg) +! merra2 data is top down +! for GFS, iflip 0: toa to sfc; 1: sfc to toa + DO ii = 1, ntrcaerm + vname=trim(specname(ii)) + call nf_inq_varid(ncid, vname, varid) + call nf_get_var(ncid, varid, buffx) + + 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 = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aerin(i,j,k,ii,2) = 1.d0*buffx(i,j,klev,1) + if(aerin(i,j,k,ii,2) < 0 .or. aerin(i,j,k,ii,2) > 1.) then + aerin(i,j,k,ii,2) = 1.e-15 + endif + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + + ENDDO ! ii-loop (ntracaerm) + +! close the file + call nf_close(ncid) + deallocate (buff, pres_tmp) + deallocate (buffx) + n1sv=n1 + n2sv=n2 + end if ! tx1 = (aer_time(n2) - rjday) / (aer_time(n2) - aer_time(n1)) tx2 = 1.0 - tx1 @@ -328,7 +515,7 @@ SUBROUTINE aerinterpol(me,master,nthrds,npts,IDATE,FHOUR,jindx1,jindx2, & !$OMP parallel num_threads(nthrds) default(none) & !$OMP shared(npts,ntrcaer,aerin,aer_pres,prsl) & !$OMP shared(ddx,ddy,jindx1,jindx2,iindx1,iindx2) & -!$OMP shared(aerpm,aerpres,aerout,n1,n2,lev,nthrds) & +!$OMP shared(aerpm,aerpres,aerout,lev,nthrds) & !$OMP shared(temij,temiy,temjx,ddxy) & !$OMP private(l,j,k,ii,i1,i2,j1,j2,tem) & !$OMP copyin(tx1,tx2) firstprivate(tx1,tx2) @@ -343,17 +530,17 @@ SUBROUTINE aerinterpol(me,master,nthrds,npts,IDATE,FHOUR,jindx1,jindx2, & I2 = IINDX2(J) DO ii=1,ntrcaer aerpm(j,L,ii) = & - tx1*(TEMIJ(j)*aerin(I1,J1,L,ii,n1)+DDXY(j)*aerin(I2,J2,L,ii,n1) & - +TEMIY(j)*aerin(I1,J2,L,ii,n1)+temjx(j)*aerin(I2,J1,L,ii,n1))& - +tx2*(TEMIJ(j)*aerin(I1,J1,L,ii,n2)+DDXY(j)*aerin(I2,J2,L,ii,n2) & - +TEMIY(j)*aerin(I1,J2,L,ii,n2)+temjx(j)*aerin(I2,J1,L,ii,n2)) + tx1*(TEMIJ(j)*aerin(I1,J1,L,ii,1)+DDXY(j)*aerin(I2,J2,L,ii,1) & + +TEMIY(j)*aerin(I1,J2,L,ii,1)+temjx(j)*aerin(I2,J1,L,ii,1))& + +tx2*(TEMIJ(j)*aerin(I1,J1,L,ii,2)+DDXY(j)*aerin(I2,J2,L,ii,2) & + +TEMIY(j)*aerin(I1,J2,L,ii,2)+temjx(j)*aerin(I2,J1,L,ii,2)) ENDDO aerpres(j,L) = & - tx1*(TEMIJ(j)*aer_pres(I1,J1,L,n1)+DDXY(j)*aer_pres(I2,J2,L,n1) & - +TEMIY(j)*aer_pres(I1,J2,L,n1)+temjx(j)*aer_pres(I2,J1,L,n1))& - +tx2*(TEMIJ(j)*aer_pres(I1,J1,L,n2)+DDXY(j)*aer_pres(I2,J2,L,n2) & - +TEMIY(j)*aer_pres(I1,J2,L,n2)+temjx(j)*aer_pres(I2,J1,L,n2)) + tx1*(TEMIJ(j)*aer_pres(I1,J1,L,1)+DDXY(j)*aer_pres(I2,J2,L,1) & + +TEMIY(j)*aer_pres(I1,J2,L,1)+temjx(j)*aer_pres(I2,J1,L,1))& + +tx2*(TEMIJ(j)*aer_pres(I1,J1,L,2)+DDXY(j)*aer_pres(I2,J2,L,2) & + +TEMIY(j)*aer_pres(I1,J2,L,2)+temjx(j)*aer_pres(I2,J1,L,2)) ENDDO ENDDO #ifndef __GFORTRAN__ From b3f627d5f5fdc65c7e342bd32815234b3bcb0294 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Wed, 10 Nov 2021 10:29:05 -0500 Subject: [PATCH 2/6] delete some debug stuff --- physics/aerinterp.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index a99fae7b0..3aeb73eb7 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -164,7 +164,6 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) enddo n1 = n2 - 1 if (n2 > 12) n2 = n2 -12 - write(*,*)"AAA0",n1, n2, iamin, iamax, jamin, jamax !! =================================================================== !! loop thru m01 - m12 for aer/pres array !! =================================================================== @@ -204,7 +203,6 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) DO ii = 1, ntrcaerm vname=trim(specname(ii)) call nf_inq_varid(ncid, vname, varid) - write(*,*)"AAA2",vname call nf_get_var(ncid, varid, buffx) do j = jamin, jamax From 8674543e2cbdd1bfbaad4d756b12b3077231fbfd Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Wed, 10 Nov 2021 10:33:10 -0500 Subject: [PATCH 3/6] handle a case with n2=13 --- physics/aerinterp.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 3aeb73eb7..e3adac99b 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -415,6 +415,7 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, endif enddo n1 = n2 - 1 + if (n2 > 12) n2 = n2 -12 ! need to read a new month if (n1>n1sv) then DO ii = 1, ntrcaerm From 102264d54f17a412eff13cd37b52fb2bf72d5a00 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Wed, 10 Nov 2021 10:48:25 -0500 Subject: [PATCH 4/6] add a output information for reading a new month merra2 --- physics/aerinterp.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index e3adac99b..162051421 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -417,7 +417,8 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, n1 = n2 - 1 if (n2 > 12) n2 = n2 -12 ! need to read a new month - if (n1>n1sv) then + if (n1.ne.n1sv) then + if (me == master) write(*,*)"read in a new month MERRA2", n2 DO ii = 1, ntrcaerm do j = jamin, jamax do k = 1, levsaer From 1ec58381a419f8d028667f9f0dda1e3ab67df647 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 16 Nov 2021 13:55:00 -0500 Subject: [PATCH 5/6] put a few lines of code in read_netfaer --- physics/aerinterp.F90 | 291 +++++++++++++----------------------------- 1 file changed, 86 insertions(+), 205 deletions(-) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 162051421..89c993a45 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -9,7 +9,7 @@ module aerinterp implicit none - private + private read_netfaer public :: read_aerdata, setindxaer, aerinterpol,read_aerdataf @@ -100,7 +100,6 @@ END SUBROUTINE read_aerdata SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) use machine, only: kind_phys, kind_io4, kind_io8 use aerclm_def - use netcdf !--- in/out integer, intent(in) :: me, master, iflip, idate(4) @@ -108,9 +107,7 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) integer, intent(inout) :: errflg real(kind=kind_phys), intent(in) :: fhour !--- locals - integer :: ncid, varid integer :: i, j, k, n, ii, imon, klev, n1, n2 - character :: fname*50, mn*2, vname*10 logical :: file_exist integer IDAT(8),JDAT(8) real(kind=kind_phys) RINC(5), rjday @@ -119,9 +116,6 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) integer w3kindreal,w3kindint 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)) @@ -129,9 +123,6 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) endif ! allocate local working arrays - allocate (buff(lonsaer, latsaer, levsw)) - allocate (pres_tmp(lonsaer, levsw)) - allocate (buffx(lonsaer, latsaer, levsw, 1)) !! found interpolation months IDAT = 0 IDAT(1) = IDATE(4) @@ -165,132 +156,12 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) n1 = n2 - 1 if (n2 > 12) n2 = n2 -12 !! =================================================================== -!! loop thru m01 - m12 for aer/pres array -!! =================================================================== - write(mn,'(i2.2)') n1 - fname=trim("aeroclim.m"//mn//".nc") - call nf_open(fname , nf90_NOWRITE, ncid) - -! ====> construct 3-d pressure array (Pa) - call nf_inq_varid(ncid, "DELP", varid) - call nf_get_var(ncid, varid, buff) - - 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, levsw - pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) - enddo !k-loop - enddo !i-loop (lon) - -! extract pres_tmp to fill aer_pres (in Pa) - do k = 1, levsaer - if ( iflip == 0 ) then ! data from toa to sfc - klev = k - else ! data from sfc to top - klev = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aer_pres(i,j,k,1) = 1.d0*pres_tmp(i,klev) - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - -! ====> construct 4-d aerosol array (kg/kg) -! merra2 data is top down -! for GFS, iflip 0: toa to sfc; 1: sfc to toa - DO ii = 1, ntrcaerm - vname=trim(specname(ii)) - call nf_inq_varid(ncid, vname, varid) - call nf_get_var(ncid, varid, buffx) - - 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 = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aerin(i,j,k,ii,1) = 1.d0*buffx(i,j,klev,1) - if(aerin(i,j,k,ii,1) < 0 .or. aerin(i,j,k,ii,1) > 1.) then - aerin(i,j,k,ii,1) = 1.e-15 - endif - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - - ENDDO ! ii-loop (ntracaerm) - -! close the file - call nf_close(ncid) + call read_netfaer(n1, iflip, 1) + call read_netfaer(n2, iflip, 2) !! =================================================================== - write(mn,'(i2.2)') n2 - fname=trim("aeroclim.m"//mn//".nc") - call nf_open(fname , nf90_NOWRITE, ncid) - -! ====> construct 3-d pressure array (Pa) - call nf_inq_varid(ncid, "DELP", varid) - call nf_get_var(ncid, varid, buff) - - 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, levsw - pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) - enddo !k-loop - enddo !i-loop (lon) - -! extract pres_tmp to fill aer_pres (in Pa) - do k = 1, levsaer - if ( iflip == 0 ) then ! data from toa to sfc - klev = k - else ! data from sfc to top - klev = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aer_pres(i,j,k,2) = 1.d0*pres_tmp(i,klev) - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - -! ====> construct 4-d aerosol array (kg/kg) -! merra2 data is top down -! for GFS, iflip 0: toa to sfc; 1: sfc to toa - DO ii = 1, ntrcaerm - vname=trim(specname(ii)) - call nf_inq_varid(ncid, vname, varid) - call nf_get_var(ncid, varid, buffx) - - 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 = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aerin(i,j,k,ii,2) = 1.d0*buffx(i,j,klev,1) - if(aerin(i,j,k,ii,2) < 0 .or. aerin(i,j,k,ii,2) > 1.) then - aerin(i,j,k,ii,2) = 1.e-15 - endif - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - - ENDDO ! ii-loop (ntracaerm) - -! close the file - call nf_close(ncid) - n1sv=n1 - n2sv=n2 + n1sv=n1 + n2sv=n2 !--- - deallocate (buff, pres_tmp) - deallocate (buffx) END SUBROUTINE read_aerdataf ! SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, & @@ -356,14 +227,12 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, ! use machine, only: kind_phys, kind_io4, kind_io8 use aerclm_def - use netcdf implicit none integer, intent(in) :: iflip integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii, klev real(kind=kind_phys) fhour,temj, tx1, tx2,temi, tem real(kind=kind_phys), dimension(npts) :: temij,temiy,temjx,ddxy - character :: fname*50, mn*2, vname*10 ! @@ -379,10 +248,6 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, integer jdow, jdoy, jday real(4) rinc4(5) integer w3kindreal,w3kindint - integer ncid, varid - real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff - real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx - real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp ! IDAT = 0 @@ -429,71 +294,7 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, enddo !j-loop (lat) ENDDO ! ii-loop (ntracaerm) !! =================================================================== - allocate (buff(lonsaer, latsaer, levsw)) - allocate (pres_tmp(lonsaer, levsw)) - allocate (buffx(lonsaer, latsaer, levsw, 1)) - - write(mn,'(i2.2)') n2 - fname=trim("aeroclim.m"//mn//".nc") - call nf_open(fname , nf90_NOWRITE, ncid) - -! ====> construct 3-d pressure array (Pa) - call nf_inq_varid(ncid, "DELP", varid) - call nf_get_var(ncid, varid, buff) - - 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, levsw - pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) - enddo !k-loop - enddo !i-loop (lon) - -! extract pres_tmp to fill aer_pres (in Pa) - do k = 1, levsaer - if ( iflip == 0 ) then ! data from toa to sfc - klev = k - else ! data from sfc to top - klev = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aer_pres(i,j,k,2) = 1.d0*pres_tmp(i,klev) - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - -! ====> construct 4-d aerosol array (kg/kg) -! merra2 data is top down -! for GFS, iflip 0: toa to sfc; 1: sfc to toa - DO ii = 1, ntrcaerm - vname=trim(specname(ii)) - call nf_inq_varid(ncid, vname, varid) - call nf_get_var(ncid, varid, buffx) - - 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 = ( levsw - k ) + 1 - endif - do i = iamin, iamax - aerin(i,j,k,ii,2) = 1.d0*buffx(i,j,klev,1) - if(aerin(i,j,k,ii,2) < 0 .or. aerin(i,j,k,ii,2) > 1.) then - aerin(i,j,k,ii,2) = 1.e-15 - endif - enddo !i-loop (lon) - enddo !k-loop (lev) - enddo !j-loop (lat) - - ENDDO ! ii-loop (ntracaerm) - -! close the file - call nf_close(ncid) - deallocate (buff, pres_tmp) - deallocate (buffx) + call read_netfaer(n2, iflip, 2) n1sv=n1 n2sv=n2 end if @@ -585,5 +386,85 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, RETURN END SUBROUTINE aerinterpol + subroutine read_netfaer(nf, iflip,nt) + use machine, only: kind_phys, kind_io4, kind_io8 + use aerclm_def + use netcdf + integer, intent(in) :: iflip, nf, nt + integer :: ncid, varid, i,j,k,ii,klev + character :: fname*50, mn*2, vname*10 + real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff + real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx + real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp + +!! =================================================================== + allocate (buff(lonsaer, latsaer, levsw)) + allocate (pres_tmp(lonsaer, levsw)) + allocate (buffx(lonsaer, latsaer, levsw, 1)) + + write(mn,'(i2.2)') nf + fname=trim("aeroclim.m"//mn//".nc") + call nf_open(fname , nf90_NOWRITE, ncid) + +! ====> construct 3-d pressure array (Pa) + call nf_inq_varid(ncid, "DELP", varid) + call nf_get_var(ncid, varid, buff) + + 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, levsw + pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) + enddo !k-loop + enddo !i-loop (lon) + +! extract pres_tmp to fill aer_pres (in Pa) + do k = 1, levsaer + if ( iflip == 0 ) then ! data from toa to sfc + klev = k + else ! data from sfc to top + klev = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aer_pres(i,j,k,nt) = 1.d0*pres_tmp(i,klev) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + +! ====> construct 4-d aerosol array (kg/kg) +! merra2 data is top down +! for GFS, iflip 0: toa to sfc; 1: sfc to toa + DO ii = 1, ntrcaerm + vname=trim(specname(ii)) + call nf_inq_varid(ncid, vname, varid) + call nf_get_var(ncid, varid, buffx) + + 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 = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aerin(i,j,k,ii,nt) = 1.d0*buffx(i,j,klev,1) + if(aerin(i,j,k,ii,nt) < 0 .or. aerin(i,j,k,ii,nt) > 1.) then + aerin(i,j,k,ii,nt) = 1.e-15 + endif + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + + ENDDO ! ii-loop (ntracaerm) + +! close the file + call nf_close(ncid) + deallocate (buff, pres_tmp) + deallocate (buffx) + return + END SUBROUTINE read_netfaer + end module aerinterp From 850971ef433a5eb50dbc97e0ab02f9dfc4d966dc Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 16 Nov 2021 15:40:04 -0500 Subject: [PATCH 6/6] adding a debug option for reading a new month data --- physics/aerinterp.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 89c993a45..30ff97dff 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -283,7 +283,9 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, if (n2 > 12) n2 = n2 -12 ! need to read a new month if (n1.ne.n1sv) then +#ifdef DEBUG if (me == master) write(*,*)"read in a new month MERRA2", n2 +#endif DO ii = 1, ntrcaerm do j = jamin, jamax do k = 1, levsaer