*
L. Bakels 2022: - This module contains all subroutines handling the * internal storage and processing of the meteorological * data, including computation of PV and boundary * layer parameters * - The reading of the meteo data happens in windfields_mod * - The vertical coordinate transformation is done in * verttransform_mod * *
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later !***************************************************************************** ! * ! L. Bakels 2022: - This module contains all subroutines handling the * ! internal storage and processing of the meteorological * ! data, including computation of PV and boundary * ! layer parameters * ! - The reading of the meteo data happens in windfields_mod * ! - The vertical coordinate transformation is done in * ! verttransform_mod * ! * !***************************************************************************** module getfields_mod use par_mod use com_mod use windfields_mod use verttransform_mod implicit none real,allocatable,dimension(:,:,:) :: & uuh, & ! wind components in x-direction [m/s] vvh, & ! wind components in y-direction [m/s] pvh, & ! potential vorticity wwh ! wind components in y-direction [m/s] real,allocatable,dimension(:,:,:,:) :: & ! Same for nexted grids uuhn, & ! vvhn, & ! pvhn, & ! wwhn, & ! pwater ! RLT added partial pressure water vapor real,allocatable,dimension(:,:,:) :: & ! For calcpv ppml, & ! ppmk ! real,allocatable,dimension(:) :: & ! For calcpar ttlev, & ! qvlev, & ! ulev, & ! vlev, & ! zlev ! private :: obukhov,richardson,scalev,calcpar,calcpar_nest,calcpv,calcpv_nest public :: getfields contains subroutine alloc_getfields implicit none integer :: stat allocate(uuh(0:nxmax-1,0:nymax-1,nuvzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate uuh' allocate(vvh(0:nxmax-1,0:nymax-1,nuvzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate vvh' allocate(pvh(0:nxmax-1,0:nymax-1,nuvzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate pvh' allocate(wwh(0:nxmax-1,0:nymax-1,nwzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate wwh' allocate(uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numbnests),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate uuhn' allocate(vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numbnests),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate vvhn' allocate(pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numbnests),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate pvhn' allocate(wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,numbnests),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate wwhn' allocate(pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate pwater' allocate(ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate ppml,ppmk' allocate(ttlev(nuvzmax),qvlev(nuvzmax),ulev(nuvzmax),vlev(nuvzmax),zlev(nuvzmax),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate ttlev,qvlev,ulev,vlev,zlev' end subroutine alloc_getfields subroutine dealloc_getfields implicit none deallocate(uuh,vvh,pvh,wwh,uuhn,vvhn,pvhn,wwhn,pwater) deallocate(ppml,ppmk) deallocate(ttlev,qvlev,ulev,vlev,zlev) end subroutine dealloc_getfields subroutine getfields(itime,nstop) ! i o !***************************************************************************** ! * ! This subroutine manages the 3 data fields to be kept in memory. * ! During the first time step of petterssen it has to be fulfilled that the * ! first data field must have |wftime|<itime, i.e. the absolute value of * ! wftime must be smaller than the absolute value of the current time in [s].* ! The other 2 fields are the next in time after the first one. * ! Pointers (memind) are used, because otherwise one would have to resort the* ! wind fields, which costs a lot of computing time. Here only the pointers * ! are resorted. * ! * ! Author: A. Stohl * ! * ! 29 April 1994 * ! * !***************************************************************************** ! Changes, Bernd C. Krueger, Feb. 2001: * ! Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block. * ! Function of nstop extended. * ! * ! Unified ECMWF and GFS builds * ! Marian Harustak, 12.5.2017 * ! - Added passing of metdata_format as it was needed by called routines * !***************************************************************************** ! * ! Variables: * ! lwindinterval [s] time difference between the two wind fields read in * ! indj indicates the number of the wind field to be read in * ! indmin remembers the number of wind fields already treated * ! memind(2) pointer, on which place the wind fields are stored * ! memtime(2) [s] times of the wind fields, which are kept in memory * ! itime [s] current time since start date of trajectory calcu- * ! lation * ! nstop > 0, if trajectory has to be terminated * ! nx,ny,nuvz,nwz field dimensions in x,y and z direction * ! uu(0:nxmax,0:nymax,nuvzmax,2) wind components in x-direction [m/s] * ! vv(0:nxmax,0:nymax,nuvzmax,2) wind components in y-direction [m/s] * ! ww(0:nxmax,0:nymax,nwzmax,2) wind components in z-direction [deltaeta/s]* ! tt(0:nxmax,0:nymax,nuvzmax,2) temperature [K] * ! ps(0:nxmax,0:nymax,2) surface pressure [Pa] * ! metdata_format format of metdata (ecmwf/gfs) * ! * ! Constants: * ! idiffmax maximum allowable time difference between 2 wind * ! fields * ! * !***************************************************************************** use class_gribfile_mod use wetdepo_mod implicit none integer :: indj,itime,nstop,memaux integer :: kz,ix,jy integer :: indmin = 1 ! Check, if wind fields are available for the current time step !************************************************************** nstop=0 if ((ldirect*wftime(1).gt.ldirect*itime).or. & (ldirect*wftime(numbwf).lt.ldirect*itime)) then write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.' write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.' nstop=4 return endif if ((ldirect*memtime(1).le.ldirect*itime).and. & (ldirect*memtime(2).gt.ldirect*itime)) then ! The right wind fields are already in memory -> don't do anything !***************************************************************** continue else if ((ldirect*memtime(2).le.ldirect*itime).and. & (memtime(2).ne.999999999)) then ! Current time is after 2nd wind field ! -> Resort wind field pointers, so that current time is between 1st and 2nd !*************************************************************************** memaux=memind(1) memind(1)=memind(2) memind(2)=memaux memtime(1)=memtime(2) ! Read a new wind field and store it on place memind(2) !****************************************************** do indj=indmin,numbwf-1 if (ldirect*wftime(indj+1).gt.ldirect*itime) then if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_temp = (count_clock - count_clock0)/real(count_rate) call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh) call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_readwind = s_readwind + ((count_clock - count_clock0)/real(count_rate)-s_temp) else call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh) end if call readwind_nest(indj+1,memind(2),uuhn,vvhn,wwhn) call calcpar(memind(2)) call calcpar_nest(memind(2)) if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh) else call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh) end if call verttransform_nest(memind(2),uuhn,vvhn,wwhn,pvhn) memtime(2)=wftime(indj+1) nstop = 1 exit endif end do indmin=indj if (WETBKDEP) then call writeprecip(itime,memind(1)) endif if ((DRYDEP).or.(lnetcdfout.eq.0)) then !$OMP PARALLEL PRIVATE(ix,jy,kz) !$OMP DO ! RLT calculate dry air density do kz=1,nuvz do jy=0,nymin1 do ix=0,nxmin1 pwater(ix,jy,kz,memind(1))=qv(ix,jy,kz,memind(1))*prs(ix,jy,kz,memind(1))/ & ((r_air/r_water)*(1.-qv(ix,jy,kz,memind(1)))+qv(ix,jy,kz,memind(1))) pwater(ix,jy,kz,memind(2))=qv(ix,jy,kz,memind(2))*prs(ix,jy,kz,memind(2))/ & ((r_air/r_water)*(1.-qv(ix,jy,kz,memind(2)))+qv(ix,jy,kz,memind(2))) rho_dry(ix,jy,kz,memind(1))=(prs(ix,jy,kz,memind(1))-pwater(ix,jy,kz,memind(1)))/ & (r_air*tt(ix,jy,kz,memind(1))) rho_dry(ix,jy,kz,memind(2))=(prs(ix,jy,kz,memind(2))-pwater(ix,jy,kz,memind(2)))/ & (r_air*tt(ix,jy,kz,memind(2))) end do end do end do ! pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv) ! rho_dry=(prs-pwater)/(r_air*tt) !$OMP END DO !$OMP END PARALLEL endif else ! No wind fields, which can be used, are currently in memory ! -> read both wind fields !*********************************************************** do indj=indmin,numbwf-1 if ((ldirect*wftime(indj).le.ldirect*itime).and. & (ldirect*wftime(indj+1).gt.ldirect*itime)) then memind(1)=1 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_temp = (count_clock - count_clock0)/real(count_rate) call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh) call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_readwind = s_readwind + ((count_clock - count_clock0)/real(count_rate)-s_temp) else call readwind_gfs(indj,memind(1),uuh,vvh,wwh) end if call readwind_nest(indj,memind(1),uuhn,vvhn,wwhn) call calcpar(memind(1)) call calcpar_nest(memind(1)) if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh) else call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh) end if call verttransform_nest(memind(1),uuhn,vvhn,wwhn,pvhn) memtime(1)=wftime(indj) memind(2)=2 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_temp = (count_clock - count_clock0)/real(count_rate) call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh) call SYSTEM_CLOCK(count_clock, count_rate, count_max) s_readwind = s_readwind + ((count_clock - count_clock0)/real(count_rate)-s_temp) else call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh) end if call readwind_nest(indj+1,memind(2),uuhn,vvhn,wwhn) call calcpar(memind(2)) call calcpar_nest(memind(2)) if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh) else call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh) end if call verttransform_nest(memind(2),uuhn,vvhn,wwhn,pvhn) memtime(2)=wftime(indj+1) nstop = 1 exit endif end do indmin=indj if (WETBKDEP) then call writeprecip(itime,memind(1)) endif if ((DRYDEP).or.(lnetcdfout.eq.0)) then !$OMP PARALLEL PRIVATE(ix,jy,kz) !$OMP DO ! RLT calculate dry air density do kz=1,nuvz do jy=0,nymin1 do ix=0,nxmin1 pwater(ix,jy,kz,memind(1))=qv(ix,jy,kz,memind(1))*prs(ix,jy,kz,memind(1))/ & ((r_air/r_water)*(1.-qv(ix,jy,kz,memind(1)))+qv(ix,jy,kz,memind(1))) pwater(ix,jy,kz,memind(2))=qv(ix,jy,kz,memind(2))*prs(ix,jy,kz,memind(2))/ & ((r_air/r_water)*(1.-qv(ix,jy,kz,memind(2)))+qv(ix,jy,kz,memind(2))) rho_dry(ix,jy,kz,memind(1))=(prs(ix,jy,kz,memind(1))-pwater(ix,jy,kz,memind(1)))/ & (r_air*tt(ix,jy,kz,memind(1))) rho_dry(ix,jy,kz,memind(2))=(prs(ix,jy,kz,memind(2))-pwater(ix,jy,kz,memind(2)))/ & (r_air*tt(ix,jy,kz,memind(2))) end do end do end do ! pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv) ! rho_dry=(prs-pwater)/(r_air*tt) !$OMP END DO !$OMP END PARALLEL endif end if lwindinterv=abs(memtime(2)-memtime(1)) if (lwindinterv.gt.idiffmax) nstop=3 end subroutine getfields subroutine calcpv(n) ! i i i o !***************************************************************************** ! * ! Calculation of potential vorticity on 3-d grid. * ! * ! Author: P. James * ! 3 February 2000 * ! * ! Adaptation to FLEXPART, A. Stohl, 1 May 2000 * ! * !***************************************************************************** ! * ! Variables: * ! n temporal index for meteorological fields (1 to 2) * ! * ! Constants: * ! * !***************************************************************************** implicit none integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr integer :: nlck real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt real :: pvavr real :: thup,thdn real,parameter :: eps=1.e-5, p0=101325 ! Set number of levels to check for adjacent theta nlck=nuvz/3 ! ! Loop over entire grid !********************** do kl=1,nuvz do jy=0,nymin1 do ix=0,nxmin1 ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n) enddo enddo enddo ! ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa !$OMP PARALLEL PRIVATE(jy,ix,kl,phi,f,tanphi,cosphi,jyvp,jyvm,jumpy,juy, & !$OMP ixvp,ixvm,jumpx,ivrp,ivrm,jux,theta,klvrp,klvrm,klpt,thetap,thetam,dthetadp, & !$OMP ii,i,ivr,kdn,kch,kup,thdn,thup,dt1,dt2,dt,vx,k,dvdx, & !$OMP jj,j,uy,dudy) !$OMP DO SCHEDULE(dynamic,1) do jy=0,nymin1 if (sglobal.and.jy.eq.0) cycle if (nglobal.and.jy.eq.nymin1) cycle ! do kl=1,nuvz ! ppml(0:nxmin1,jy,kl)=akz(kl)+bkz(kl)*ps(0:nxmin1,jy,1,n) ! ppmk(0:nxmin1,jy,kl)=(100000./ppml(0:nxmin1,jy,kl))**kappa ! end do phi = (ylat0 + jy * dy) * pi / 180. f = 0.00014585 * sin(phi) tanphi = tan(phi) cosphi = cos(phi) ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat) jyvp=jy+1 jyvm=jy-1 if (jy.eq.0) jyvm=0 if (jy.eq.nymin1) jyvp=nymin1 ! Define absolute gap length jumpy=2 if (jy.eq.0.or.jy.eq.nymin1) jumpy=1 if (sglobal.and.jy.eq.1) then jyvm=1 jumpy=1 end if if (nglobal.and.jy.eq.ny-2) then jyvp=ny-2 jumpy=1 end if juy=jumpy ! do ix=0,nxmin1 ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long) ixvp=ix+1 ixvm=ix-1 jumpx=2 if (xglobal) then ivrp=ixvp ivrm=ixvm if (ixvm.lt.0) ivrm=ixvm+nxmin1 if (ixvp.ge.nx) ivrp=ixvp-nx+1 else if (ix.eq.0) ixvm=0 if (ix.eq.nxmin1) ixvp=nxmin1 ivrp=ixvp ivrm=ixvm ! Define absolute gap length if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1 end if jux=jumpx ! ! Loop over the vertical !*********************** do kl=1,nuvz theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl) klvrp=kl+1 klvrm=kl-1 klpt=kl ! If top or bottom level, dthetadp is evaluated between the current ! level and the level inside, otherwise between level+1 and level-1 ! if (klvrp.gt.nuvz) klvrp=nuvz if (klvrm.lt.1) klvrm=1 thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp) thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm) dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm)) ! Compute vertical position at pot. temperature surface on subgrid ! and the wind at that position !***************************************************************** ! a) in x direction ii=0 x_loop: do i=ixvm,ixvp,jumpx ivr=i if (xglobal) then if (i.lt.0) ivr=ivr+nxmin1 if (i.ge.nx) ivr=ivr-nx+1 end if ii=ii+1 ! Search adjacent levels for current theta value ! Spiral out from current level for efficiency kup=klpt-1 kdn=klpt kch=0 x_lev_loop: do while (kch.lt.nlck) ! Upward branch kup=kup+1 if (kup.lt.nuvz) then kch=kch+1 k=kup thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k) thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt cycle x_loop endif endif ! Downward branch kdn=kdn-1 if (kdn.ge.1) then kch=kch+1 k=kdn thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k) thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt cycle x_loop endif endif end do x_lev_loop ! This section used when no values were found ! Must use vv at current level and long. jux becomes smaller by 1 vx(ii)=vvh(ix,jy,kl) jux=jux-1 ! Otherwise OK end do x_loop if (jux.gt.0) then dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.) else dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl) dvdx=dvdx/real(jumpx)/(dx*pi/180.) ! Only happens if no equivalent theta value ! can be found on either side, hence must use values ! from either side, same pressure level. end if ! b) in y direction jj=0 y_loop: do j=jyvm,jyvp,jumpy jj=jj+1 ! Search adjacent levels for current theta value ! Spiral out from current level for efficiency kup=klpt-1 kdn=klpt kch=0 y_lev_loop: do while (kch.lt.nlck) ! Upward branch kup=kup+1 if (kup.lt.nuvz) then kch=kch+1 k=kup thdn=tth(ix,j,k,n)*ppmk(ix,j,k) thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt cycle y_loop endif endif ! Downward branch kdn=kdn-1 if (kdn.ge.1) then kch=kch+1 k=kdn thdn=tth(ix,j,k,n)*ppmk(ix,j,k) thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt cycle y_loop endif endif end do y_lev_loop ! This section used when no values were found ! Must use uu at current level and lat. juy becomes smaller by 1 uy(jj)=uuh(ix,jy,kl) juy=juy-1 ! Otherwise OK end do y_loop if (juy.gt.0) then dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.) else dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl) dudy=dudy/real(jumpy)/(dy*pi/180.) end if ! pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy & +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81 ! ! Resest jux and juy jux=jumpx juy=jumpy end do end do end do !$OMP END DO !$OMP END PARALLEL ! ! Fill in missing PV values on poles, if present ! Use mean PV of surrounding latitude ring ! if (sglobal) then do kl=1,nuvz pvavr=0. do ix=0,nxmin1 pvavr=pvavr+pvh(ix,1,kl) end do pvavr=pvavr/real(nx) jy=0 do ix=0,nxmin1 pvh(ix,jy,kl)=pvavr end do end do end if if (nglobal) then do kl=1,nuvz pvavr=0. do ix=0,nxmin1 pvavr=pvavr+pvh(ix,ny-2,kl) end do pvavr=pvavr/real(nx) jy=nymin1 do ix=0,nxmin1 pvh(ix,jy,kl)=pvavr end do end do end if end subroutine calcpv subroutine calcpv_nest(l,n) ! i i i i o !***************************************************************************** ! * ! Calculation of potential vorticity on 3-d nested grid * ! * ! Author: P. James * ! 22 February 2000 * ! * !***************************************************************************** ! * ! Variables: * ! n temporal index for meteorological fields (1 to 2) * ! l index of current nest * ! * ! Constants: * ! * !***************************************************************************** implicit none integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr integer :: nlck,l real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt real :: thup,thdn real,parameter :: eps=1.e-5,p0=101325 ! Set number of levels to check for adjacent theta nlck=nuvz/3 ! ! Loop over entire grid !********************** do kl=1,nuvz do jy=0,nyn(l)-1 do ix=0,nxn(l)-1 ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l) enddo enddo enddo ! ppmk=(100000./ppml)**kappa ppmk(0:nxn(l)-1,0:nyn(l)-1,1:nuvz)=(100000./ppml(0:nxn(l)-1,0:nyn(l)-1,1:nuvz))**kappa do jy=0,nyn(l)-1 phi = (ylat0n(l) + jy * dyn(l)) * pi / 180. f = 0.00014585 * sin(phi) tanphi = tan(phi) cosphi = cos(phi) ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat) jyvp=jy+1 jyvm=jy-1 if (jy.eq.0) jyvm=0 if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1 ! Define absolute gap length jumpy=2 if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1 juy=jumpy ! do ix=0,nxn(l)-1 ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long) ixvp=ix+1 ixvm=ix-1 jumpx=2 if (ix.eq.0) ixvm=0 if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1 ivrp=ixvp ivrm=ixvm ! Define absolute gap length if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1 jux=jumpx ! ! Loop over the vertical !*********************** do kl=1,nuvz theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl) klvrp=kl+1 klvrm=kl-1 klpt=kl ! If top or bottom level, dthetadp is evaluated between the current ! level and the level inside, otherwise between level+1 and level-1 ! if (klvrp.gt.nuvz) klvrp=nuvz if (klvrm.lt.1) klvrm=1 thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp) thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm) dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm)) ! Compute vertical position at pot. temperature surface on subgrid ! and the wind at that position !***************************************************************** ! a) in x direction ii=0 x_loop: do i=ixvm,ixvp,jumpx ivr=i ii=ii+1 ! Search adjacent levels for current theta value ! Spiral out from current level for efficiency kup=klpt-1 kdn=klpt kch=0 x_lev_loop: do while (kch.lt.nlck) ! Upward branch kup=kup+1 if (kup.lt.nuvz) then kch=kch+1 k=kup thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k) thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt cycle x_loop endif endif ! Downward branch kdn=kdn-1 if (kdn.ge.1) then kch=kch+1 k=kdn thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k) thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt cycle x_loop endif endif end do x_lev_loop ! This section used when no values were found ! Must use vv at current level and long. jux becomes smaller by 1 vx(ii)=vvhn(ix,jy,kl,l) jux=jux-1 ! Otherwise OK end do x_loop if (jux.gt.0) then dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.) else dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l) dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.) ! Only happens if no equivalent theta value ! can be found on either side, hence must use values ! from either side, same pressure level. end if ! b) in y direction jj=0 y_loop: do j=jyvm,jyvp,jumpy jj=jj+1 ! Search adjacent levels for current theta value ! Spiral out from current level for efficiency kup=klpt-1 kdn=klpt kch=0 y_lev_loop: do while (kch.lt.nlck) ! Upward branch kup=kup+1 if (kup.lt.nuvz) then kch=kch+1 k=kup thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k) thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt cycle y_loop endif endif ! Downward branch kdn=kdn-1 if (kdn.ge.1) then kch=kch+1 k=kdn thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k) thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1) if (((thdn.ge.theta).and.(thup.le.theta)).or. & ((thdn.le.theta).and.(thup.ge.theta))) then dt1=abs(theta-thdn) dt2=abs(theta-thup) dt=dt1+dt2 if (dt.lt.eps) then ! Avoid division by zero error dt1=0.5 ! G.W., 10.4.1996 dt2=0.5 dt=1.0 endif uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt cycle y_loop endif endif end do y_lev_loop ! This section used when no values were found ! Must use uu at current level and lat. juy becomes smaller by 1 uy(jj)=uuhn(ix,jy,kl,l) juy=juy-1 ! Otherwise OK end do y_loop if (juy.gt.0) then dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.) else dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l) dudy=dudy/real(jumpy)/(dyn(l)*pi/180.) end if pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy & +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81 ! Resest jux and juy jux=jumpx juy=jumpy end do end do end do end subroutine calcpv_nest subroutine calcpar(n) ! i i i o !***************************************************************************** ! * ! Computation of several boundary layer parameters needed for the * ! dispersion calculation and calculation of dry deposition velocities. * ! All parameters are calculated over the entire grid. * ! * ! Author: A. Stohl * ! * ! 21 May 1995 * ! * !***************************************************************************** ! Petra Seibert, Feb 2000: * ! convection scheme: * ! new variables in call to richardson * ! * ! CHANGE 17/11/2005 Caroline Forster NCEP GFS version * ! * ! Changes, Bernd C. Krueger, Feb. 2001: * ! Variables tth and qvh (on eta coordinates) in common block * ! * ! Unified ECMWF and GFS builds * ! Marian Harustak, 12.5.2017 * ! - Merged calcpar and calcpar_gfs into one routine using if-then * ! for meteo-type dependent code * !***************************************************************************** ! Changes Anne Tipka June 2023: * ! sum up precipitation fields over number of available fields in a single * ! time interval (newWetDepoScheme) * !***************************************************************************** ! * ! Variables: * ! n temporal index for meteorological fields (1 to 3) * ! metdata_format format of metdata (ecmwf/gfs) * ! * ! Constants: * ! * ! * ! Functions: * ! scalev computation of ustar * ! obukhov computatio of Obukhov length * ! * !***************************************************************************** use class_gribfile_mod use drydepo_mod use qvsat_mod implicit none integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start,ierr,stat real :: ol,hmixplus real :: rh,subsceff,ylat real :: altmin,tvold,pold,zold,pint,tv,hmixdummy,akzdummy real,allocatable,dimension(:) :: vd real,parameter :: const=r_air/ga !write(*,*) 'in calcpar writting snowheight' !*********************************** ! for test: write out snow depths ! open(4,file='slandusetest',form='formatted') ! do 5 ix=0,nxmin1 !5 write (4,*) (sd(ix,jy,1,n),jy=0,nymin1) ! close(4) ! Loop over entire grid !********************** !$OMP PARALLEL PRIVATE(jy,ix,ulev,vlev,ttlev,qvlev,llev,ylat,ol,i,hmixplus, & !$OMP subsceff,vd,kz,lz,zlev,rh,kzmin,pold,zold,tvold,pint,tv,loop_start,ierr, & !$OMP altmin) allocate( vd(maxspec),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate vd inside of OMP loop' !$OMP DO do jy=0,nymin1 ! Set minimum height for tropopause !********************************** ylat=ylat0+real(jy)*dy if ((ylat.ge.-20.).and.(ylat.le.20.)) then altmin = 5000. else if ((ylat.gt.20.).and.(ylat.lt.40.)) then altmin=2500.+(40.-ylat)*125. else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then altmin=2500.+(40.+ylat)*125. else altmin=2500. endif endif do ix=0,nxmin1 ! 1) Calculation of friction velocity !************************************ ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), & td2(ix,jy,1,n),sfcstress(ix,jy,1,n)) if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8 ! 2) Calculation of inverse Obukhov length scale !*********************************************** if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then ! NCEP version: find first level above ground llev = 0 do i=1,nuvz if (ps(ix,jy,1,n).lt.akz(i)) llev=i end do llev = llev+1 if (llev.gt.nuvz) llev = nuvz-1 ! NCEP version ! calculate inverse Obukhov length scale with tth(llev) ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), & tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n), & akm,bkm,akz(llev)) else llev=0 ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), & tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akzdummy) end if if (ol.ne.0.) then oli(ix,jy,1,n)=1./ol else oli(ix,jy,1,n)=99999. endif ! 3) Calculation of convective velocity scale and mixing height !************************************************************** do i=1,nuvz ulev(i)=uuh(ix,jy,i) vlev(i)=vvh(ix,jy,i) ttlev(i)=tth(ix,jy,i,n) qvlev(i)=qvh(ix,jy,i,n) end do if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, & ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), & td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,ierr) else call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, & ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), & td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,ierr) end if if (ierr.lt.0) then write(*,9500) 'failure', ix, jy error stop 'calcpar: richardson computation failed' endif 9500 format( 'calcpar - richardson ', a, ' - ix,jy=', 2i5 ) if(lsubgrid.eq.1) then subsceff=min(excessoro(ix,jy),hmixplus) else subsceff=0.0 endif ! ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY ! hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height ! 4) Calculation of dry deposition velocities !******************************************** if (DRYDEP) then ! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on ! windspeed z0_drydep(ix,jy)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga ! Calculate relative humidity at surface !*************************************** rh=ew(td2(ix,jy,1,n),ps(ix,jy,1,n))/ew(tt2(ix,jy,1,n),ps(ix,jy,1,n)) call getvdep(n,ix,jy,ustar(ix,jy,1,n), & tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), & ssr(ix,jy,1,n),rh,sum(lsprec(ix,jy,1,:,n))+sum(convprec(ix,jy,1,:,n)), & sd(ix,jy,1,n),vd) do i=1,nspec vdep(ix,jy,i,n)=vd(i) end do endif !****************************************************** ! Calculate height of thermal tropopause (Hoinka, 1997) !****************************************************** ! 1) Calculate altitudes of model levels !*************************************** tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n),ps(ix,jy,1,n))/ & ps(ix,jy,1,n)) pold=ps(ix,jy,1,n) zold=0. if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then loop_start=llev else loop_start=2 end if do kz=loop_start,nuvz pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n)) if (abs(tv-tvold).gt.0.2) then zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold) else zlev(kz)=zold+const*log(pold/pint)*tv endif tvold=tv pold=pint zold=zlev(kz) end do ! 2) Define a minimum level kzmin, from which upward the tropopause is ! searched for. This is to avoid inversions in the lower troposphere ! to be identified as the tropopause !************************************************************************ if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then !LB, The CTM version has 2 (as bugfix), so I changed it 2 from 1 to try out loop_start=2 else loop_start=llev end if kzmin=nuvz do kz=loop_start,nuvz if (zlev(kz).ge.altmin) then kzmin=kz exit endif end do ! 3) Search for first stable layer above minimum height that fulfills the ! thermal tropopause criterion !************************************************************************ outer: do kz=kzmin,nuvz inner: do lz=kz+1,nuvz if ((zlev(lz)-zlev(kz)).gt.2000.) then if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ & (zlev(lz)-zlev(kz))).lt.0.002) then tropopause(ix,jy,1,n)=zlev(kz) exit outer endif exit inner endif end do inner end do outer end do end do !$OMP END DO deallocate(vd) !$OMP END PARALLEL ! Calculation of potential vorticity on 3-d grid !*********************************************** call calcpv(n) end subroutine calcpar subroutine calcpar_nest(n) ! i i i o !***************************************************************************** ! * ! Computation of several boundary layer parameters needed for the * ! dispersion calculation and calculation of dry deposition velocities. * ! All parameters are calculated over the entire grid. * ! This routine is similar to calcpar, but is used for the nested grids. * ! * ! Author: A. Stohl * ! * ! 8 February 1999 * ! * !***************************************************************************** ! Petra Seibert, Feb 2000: * ! convection scheme: * ! new variables in call to richardson * ! * !***************************************************************************** ! Changes, Bernd C. Krueger, Feb. 2001: * ! Variables tth and qvh (on eta coordinates) in common block * ! * ! Unified ECMWF and GFS builds * ! Marian Harustak, 12.5.2017 * !***************************************************************************** ! Changes Anne Tipka June 2023: * ! sum up precipitation fields over number of available fields in a single * ! time interval (newWetDepoScheme) * !***************************************************************************** ! * ! Variables: * ! n temporal index for meteorological fields (1 to 3) * ! metdata_format format of metdata (ecmwf/gfs) * ! * ! Constants: * ! * ! * ! Functions: * ! scalev computation of ustar * ! obukhov computatio of Obukhov length * ! * !***************************************************************************** use drydepo_mod use qvsat_mod implicit none integer :: n,ix,jy,i,l,kz,lz,kzmin,ierr,stat real :: ol,hmixplus,dummyakzllev real :: rh,subsceff,ylat real :: altmin,tvold,pold,zold,pint,tv real,allocatable,dimension(:) :: vd real,parameter :: const=r_air/ga ! Loop over all nests !******************** do l=1,numbnests ! Loop over entire grid !********************** !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i,ix,jy,kz,lz,kzmin,tvold,pold,zold,zlev,tv,pint, & !$OMP rh,ierr,subsceff,ulev,vlev,ttlev,qvlev,ol,altmin,ylat,hmixplus, & !$OMP dummyakzllev,stat,vd ) allocate( vd(maxspec),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate vd inside of OMP loop' !$OMP DO do jy=0,nyn(l)-1 ! Set minimum height for tropopause !********************************** ylat=ylat0n(l)+real(jy)*dyn(l) if ((ylat.ge.-20.).and.(ylat.le.20.)) then altmin = 5000. else if ((ylat.gt.20.).and.(ylat.lt.40.)) then altmin=2500.+(40.-ylat)*125. else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then altmin=2500.+(40.+ylat)*125. else altmin=2500. endif endif do ix=0,nxn(l)-1 ! 1) Calculation of friction velocity !************************************ ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), & td2n(ix,jy,1,n,l),sfcstressn(ix,jy,1,n,l)) if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8 ! 2) Calculation of inverse Obukhov length scale !*********************************************** ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), & td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), & sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev) if (ol.ne.0.) then olin(ix,jy,1,n,l)=1./ol else olin(ix,jy,1,n,l)=99999. endif ! 3) Calculation of convective velocity scale and mixing height !************************************************************** do i=1,nuvz ulev(i)=uuhn(ix,jy,i,l) vlev(i)=vvhn(ix,jy,i,l) ttlev(i)=tthn(ix,jy,i,n,l) qvlev(i)=qvhn(ix,jy,i,n,l) end do call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, & qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), & tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), & wstarn(ix,jy,1,n,l),hmixplus,ierr) if (ierr.lt.0) then write(*,9500) 'failure', ix, jy, l error stop 'calcpar_nest: richardson computation failed' endif 9500 format( 'calcparn - richardson ', a, ' - ix,jy=', 2i5 ) if(lsubgrid.eq.1) then subsceff=min(excessoron(ix,jy,l),hmixplus) else subsceff=0.0 endif ! ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY ! hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height ! 4) Calculation of dry deposition velocities !******************************************** if (DRYDEP) then ! z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga ! z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga z0_drydepn(ix,jy,l)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga ! Calculate relative humidity at surface !*************************************** rh=ew(td2n(ix,jy,1,n,l),psn(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l)) call getvdep_nest(n,ix,jy,ustarn(ix,jy,1,n,l), & tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), & ssrn(ix,jy,1,n,l),rh,sum(lsprecn(ix,jy,1,:,n,l))+ & sum(convprecn(ix,jy,1,:,n,l)),sdn(ix,jy,1,n,l),vd,l) do i=1,nspec vdepn(ix,jy,i,n,l)=vd(i) end do endif !****************************************************** ! Calculate height of thermal tropopause (Hoinka, 1997) !****************************************************** ! 1) Calculate altitudes of ECMWF model levels !********************************************* tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l),psn(ix,jy,1,n,l))/ & psn(ix,jy,1,n,l)) pold=psn(ix,jy,1,n,l) zold=0. do kz=2,nuvz pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) ! pressure on model layers tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l)) if (abs(tv-tvold).gt.0.2) then zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold) else zlev(kz)=zold+const*log(pold/pint)*tv endif tvold=tv pold=pint zold=zlev(kz) end do ! 2) Define a minimum level kzmin, from which upward the tropopause is ! searched for. This is to avoid inversions in the lower troposphere ! to be identified as the tropopause !************************************************************************ kzmin=1 do kz=1,nuvz if (zlev(kz).ge.altmin) then kzmin=kz exit endif end do ! 3) Search for first stable layer above minimum height that fulfills the ! thermal tropopause criterion !************************************************************************ kzloop : do kz=kzmin,nuvz lzloop : do lz=kz+1,nuvz if ((zlev(lz)-zlev(kz)).gt.2000.) then if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ & (zlev(lz)-zlev(kz))).lt.0.002) then tropopausen(ix,jy,1,n,l)=zlev(kz) exit kzloop endif exit lzloop endif end do lzloop end do kzloop end do end do !$OMP END DO deallocate(vd) !$OMP END PARALLEL ! Calculation of potential vorticity on 3-d grid !*********************************************** call calcpv_nest(l,n) end do end subroutine calcpar_nest real function obukhov(ps,tsfc,tdsfc,tlev,ustar,hf,akm,bkm,plev) !******************************************************************** ! * ! Author: G. WOTAWA * ! Date: 1994-06-27 * ! * ! This program calculates Obukhov scale height from surface * ! meteorological data and sensible heat flux. * ! * !******************************************************************** ! * ! Update: A. Stohl, 2000-09-25, avoid division by zero by * ! setting ustar to minimum value * ! CHANGE: 17/11/2005 Caroline Forster NCEP GFS version * ! * ! Unified ECMWF and GFS builds * ! Marian Harustak, 12.5.2017 * ! - Merged obukhov and obukhov_gfs into one routine using * ! if-then for meteo-type dependent code * ! * !******************************************************************** ! * ! INPUT: * ! * ! ps surface pressure [Pa] * ! tsfc surface temperature [K] * ! tdsfc surface dew point [K] * ! tlev temperature first model level [K] * ! ustar scale velocity [m/s] * ! hf surface sensible heat flux [W/m2] * ! akm ECMWF vertical discretization parameter * ! bkm ECMWF vertical discretization parameter * ! plev * ! metdata_format format of metdata (ecmwf/gfs) * ! * !******************************************************************** use class_gribfile_mod use qvsat_mod implicit none real,dimension(:) :: akm,bkm real :: ps,tsfc,tdsfc,tlev,ustar,hf,e,tv,rhoa,plev real :: ak1,bk1,theta,thetastar e=ew(tdsfc,ps) ! vapor pressure tv=tsfc*(1.+0.378*e/ps) ! virtual temperature rhoa=ps/(r_air*tv) ! air density if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then ak1=(akm(1)+akm(2))*0.5 bk1=(bkm(1)+bkm(2))*0.5 plev=ak1+bk1*ps ! Pressure level 1 end if theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature if (ustar.le.0.) ustar=1.e-8 thetastar=hf/(rhoa*cpa*ustar) ! scale temperature if(abs(thetastar).gt.1.e-10) then obukhov=theta*ustar**2/(karman*ga*thetastar) else obukhov=9999 ! zero heat flux endif if (obukhov.gt. 9999.) obukhov= 9999. if (obukhov.lt.-9999.) obukhov=-9999. end function obukhov subroutine richardson(psfc,ust,ttlev,qvlev,ulev,vlev,nuvz, & akz,bkz,hf,tt2,td2,h,wst,hmixplus,ierr) ! i i i i i i i ! i i i i i o o o !**************************************************************************** ! * ! Calculation of mixing height based on the critical Richardson number. * ! Calculation of convective time scale. * ! For unstable conditions, one iteration is performed. An excess * ! temperature (dependent on hf and wst) is calculated, added to the * ! temperature at the lowest model level. Then the procedure is repeated.* ! * ! Author: A. Stohl * ! * ! 22 August 1996 * ! * ! Literature: * ! Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts * ! of alternative boundary-layer height formulations. Boundary-Layer * ! Meteor. 81, 245-269. * ! * !**************************************************************************** ! * ! Update: 1999-02-01 by G. Wotawa * ! * ! Two meter level (temperature, humidity) is taken as reference level * ! instead of first model level. * ! New input variables tt2, td2 introduced. * ! * ! CHANGE: 17/11/2005 Caroline Forster NCEP GFS version * ! * ! Unified ECMWF and GFS builds * ! Marian Harustak, 12.5.2017 * ! - Merged richardson and richardson_gfs into one routine using * ! if-then for meteo-type dependent code * ! * !**************************************************************************** ! * ! Variables: * ! h mixing height [m] * ! hf sensible heat flux * ! psfc surface pressure at point (xt,yt) [Pa] * ! tv virtual temperature * ! wst convective velocity scale * ! metdata_format format of metdata (ecmwf/gfs) * ! * ! Constants: * ! ric critical Richardson number * ! * !**************************************************************************** use class_gribfile_mod use qvsat_mod implicit none integer,intent(out) :: & ierr ! Returns error when no richardson number can be found real, intent(out) :: & h, & ! mixing height [m] wst, & ! convective velocity scale hmixplus ! integer,intent(in) :: & nuvz ! Upper vertical level real,intent(in) :: & psfc, & ! surface pressure at point (xt,yt) [Pa] ust, & ! Scale velocity hf, & ! Surface sensible heat flux tt2,td2 ! Temperature real,intent(in),dimension(:) :: & ttlev, & qvlev, & ulev, & vlev, & akz,bkz integer :: & i,k,iter,llev,loop_start,kcheck ! Loop variables real :: & tv,tvold, & ! Virtual temperature zref,z,zold,zl,zl1,zl2, & ! Heights pint,pold, & ! Pressures theta,thetaold,thetaref,thetal, & ! Potential temperature theta1,theta2,thetam, & ri, & ! Richardson number per level ril, & ! Richardson number sub level excess, & ! ul,vl, & ! Velocities sub level wspeed, & ! Wind speed at z=hmix bvfsq, & ! Brunt-Vaisala frequency bvf, & ! square root of bvfsq rh,rhold,rhl real,parameter :: const=r_air/ga, ric=0.25, b=100., bs=8.5 integer,parameter :: itmax=3 excess=0.0 if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then ! NCEP version: find first model level above ground !************************************************** llev = 0 do i=1,nuvz if (psfc.lt.akz(i)) llev=i end do llev = llev+1 ! sec llev should not be 1! if (llev.eq.1) llev = 2 if (llev.gt.nuvz) llev = nuvz-1 ! NCEP version end if ! Compute virtual temperature and virtual potential temperature at ! reference level (2 m) !***************************************************************** do iter=1,itmax,1 pold=psfc tvold=tt2*(1.+0.378*ew(td2,psfc)/psfc) zold=2.0 zref=zold rhold=ew(td2,psfc)/ew(tt2,psfc) thetaref=tvold*(100000./pold)**(r_air/cpa)+excess thetaold=thetaref ! Integrate z up to one level above zt !************************************* if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then loop_start=llev else loop_start=2 end if kcheck=loop_start do k=loop_start,nuvz kcheck=k pint=akz(k)+bkz(k)*psfc ! pressure on model layers tv=ttlev(k)*(1.+0.608*qvlev(k)) if (abs(tv-tvold).gt.0.2) then z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold) else z=zold+const*log(pold/pint)*tv endif theta=tv*(100000./pint)**(r_air/cpa) ! PS rh = qvlev(k) / f_qvsat( pint, ttlev(k) ) ! Calculate Richardson number at each level !**************************************** ri=ga/thetaref*(theta-thetaref)*(z-zref)/ & max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1) ! addition of second condition: MH should not be placed in an ! unstable layer (PS / Feb 2000) if (ri.gt.ric .and. thetaold.lt.theta) exit tvold=tv pold=pint rhold=rh thetaold=theta zold=z end do ! Check opied from FLEXPART-WRF, 2022 LB if (kcheck.ge.nuvz) then write(*,*) 'richardson not working, no stable layer -- k = nuvz' ierr = -10 goto 7000 endif !k=min(k,nuvz) ! ESO: make sure k <= nuvz (ticket #139) !MD change to work without goto ! Determine Richardson number between the critical levels !******************************************************** zl1=zold theta1=thetaold do i=1,20 zl=zold+real(i)/20.*(z-zold) ul=ulev(kcheck-1)+real(i)/20.*(ulev(kcheck)-ulev(kcheck-1)) vl=vlev(kcheck-1)+real(i)/20.*(vlev(kcheck)-vlev(kcheck-1)) thetal=thetaold+real(i)/20.*(theta-thetaold) rhl=rhold+real(i)/20.*(rh-rhold) ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ & max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1) zl2=zl theta2=thetal if (ril.gt.ric) exit if (i.eq.20) then write(*,*) 'WARNING: NO RICHARDSON NUMBER GREATER THAN 0.25 FOUND', kcheck,nuvz,ril,ri exit endif zl1=zl theta1=thetal if (i.eq.20) error stop 'RICHARDSON: NO RICHARDSON NUMBER GREATER THAN 0.25 FOUND' end do h=zl thetam=0.5*(theta1+theta2) wspeed=sqrt(ul**2+vl**2) ! Wind speed at z=hmix bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1) ! Brunt-Vaisala frequency ! at z=hmix ! Under stable conditions, limit the maximum effect of the subgrid-scale topography ! by the maximum lifting possible from the available kinetic energy !***************************************************************************** if(bvfsq.le.0.) then hmixplus=9999. else bvf=sqrt(bvfsq) hmixplus=wspeed/bvf*convke endif ! Calculate convective velocity scale !************************************ if (hf.lt.0.) then wst=(-h*ga/thetaref*hf/cpa)**0.333 excess=-bs*hf/cpa/wst else wst=0. exit endif end do ierr = 0 return ! Fatal error -- print the inputs 7000 continue write(*,'(a )') 'nuvz' write(*,'(i5 )') nuvz write(*,'(a )') 'psfc,ust,hf,tt2,td2,h,wst,hmixplus' write(*,'(1p,4e18.10)') psfc,ust,hf,tt2,td2,h,wst,hmixplus return end subroutine richardson real function scalev(ps,t,td,stress) !******************************************************************** ! * ! Author: G. WOTAWA * ! Date: 1994-06-27 * ! Update: 1996-05-21 A. Stohl * ! * !******************************************************************** ! * ! This Programm calculates scale velocity ustar from surface * ! stress and air density. * ! * !******************************************************************** ! * ! INPUT: * ! * ! ps surface pressure [Pa] * ! t surface temperature [K] * ! td surface dew point [K] * ! stress surface stress [N/m2] * ! * !******************************************************************** use qvsat_mod implicit none real :: ps,t,td,e,tv,rhoa,stress e=ew(td,ps) ! vapor pressure tv=t*(1.+0.378*e/ps) ! virtual temperature rhoa=ps/(r_air*tv) ! air density scalev=sqrt(abs(stress)/rhoa) end function scalev end module getfields_mod