diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index 94fc5e36b..57d253083 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -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 @@ -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 @@ -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) & @@ -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, & @@ -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 @@ -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 diff --git a/physics/aerclm_def.F b/physics/aerclm_def.F index 84852a1de..426881fe4 100644 --- a/physics/aerclm_def.F +++ b/physics/aerclm_def.F @@ -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) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index e7cd6ca20..bed73c5be 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -11,7 +11,7 @@ module aerinterp private - public :: read_aerdata, setindxaer, aerinterpol + public :: read_aerdata, setindxaer, aerinterpol, read_aerdataf contains @@ -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 @@ -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 !! =================================================================== @@ -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) @@ -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) @@ -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) @@ -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) @@ -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 @@ -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 diff --git a/physics/radiation_aerosols.f b/physics/radiation_aerosols.f index f732c37ef..e1e66b0d9 100644 --- a/physics/radiation_aerosols.f +++ b/physics/radiation_aerosols.f @@ -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. @@ -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) diff --git a/physics/sfc_sice.f b/physics/sfc_sice.f index ab67f849e..081bbf48e 100644 --- a/physics/sfc_sice.f +++ b/physics/sfc_sice.f @@ -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) @@ -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