!################################################################# !BOP ! ! !MODULE: DIFFUSION ! ! !DESCRIPTION: ! ! Diffusion routines. ! ! ! !REVISION HISTORY: ----------------------------- ! ! ????-?? Bram Bregman. ! Adjusted for TM5 by ! ! 2005-04 ?? ! Adjoint version ! ! 2010-08 Arjo Segers ! Debugged computation of u* over sea: ! abs. surface stress should be devided by air density. ! Set minimum value for u* to trap division by zero (MK). ! ! 2012-01 Philippe Le Sager ! Adapted for Lat-Lon MPI domain decomposition ! ! 2019-06 Andy Jacobson ! Changed DKG file format to CF-compliant netCDF. ! Also added diffusion_filename function (Sourish Basu) ! to place DKG files in a meteo file hierarchy. ! !EOP !################################################################# ! ! global input fields: ! phlb(:,:,lm+1) ! half level pressure [Pa] ! gph(:,:,lm+1) ! height of half-levels [m] !CMK note gph starts NOW at 1! ! t ! temperature [K] ! q ! specific humidity [kg/kg] ! uwind ! mean wind W-E [m/s] ! vwind ! mean wind S-N [m/s] ! sshf ! surface sensible heat flux [W / m^2] ! slhf ! surface latent heat flux [W / m^2] ! ustar (ustar_loc) ! velocity scale ! ! global output fields: ! kvh ! eddy diffusivity for heat at top [m^2/s] ! pblh ! boundary layer height[m] ! dkg ! vertical diffusion coefficient [ks s-1] ! !### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# MODULE DIFFUSION use GO, only : gol, goPr, goErr use dims, only : nregions use global_types, only : emis_data implicit none private ! methods public :: Diffusion_Init, Diffusion_Done public :: Calc_Kzz, Read_Diffusion, Write_Diffusion ! flag public :: diffusion_write ! --- const -------------------------------------- character(len=*), parameter :: mname = 'diffusion' logical :: diffusion_write = .false. character(len=2000) :: diffusion_dir = '' character(len=20) :: mlevs ! --- local ------------------------------------ type(emis_data), dimension(nregions), target :: ustar_loc, sr_mix CONTAINS ! ================================================================ subroutine Diffusion_Init( status ) use global_data, only : rcfile use MeteoData, only : Set use MeteoData, only : pu_dat, pv_dat use MeteoData, only : temper_dat, humid_dat, gph_dat use MeteoData, only : oro_dat, lsmask_dat use MeteoData, only : sr_ecm_dat, sr_ols_dat use MeteoData, only : wspd_dat use MeteoData, only : slhf_dat, sshf_dat use MeteoData, only : ewss_dat, nsss_dat use GO, only : TrcFile, Init, Done, ReadRc ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Diffusion_Init' ! --- local ------------------------------------- integer :: region type(TrcFile) :: rcF ! --- begin -------------------------------- ! loop over all (zoom) regions: do region = 1, nregions ! meteo used by diffusion call Set( pu_dat(region), status, used=.true. ) call Set( pv_dat(region), status, used=.true. ) call Set( temper_dat(region), status, used=.true. ) call Set( humid_dat(region), status, used=.true. ) call Set( gph_dat(region), status, used=.true. ) call Set( oro_dat(region), status, used=.true. ) call Set( lsmask_dat(region), status, used=.true. ) call Set( sr_ecm_dat(region), status, used=.true. ) call Set( sr_ols_dat(region), status, used=.true. ) call Set( wspd_dat(region), status, used=.true. ) call Set( sshf_dat(region), status, used=.true. ) call Set( slhf_dat(region), status, used=.true. ) call Set( ewss_dat(region), status, used=.true. ) call Set( nsss_dat(region), status, used=.true. ) end do call Init( rcF, rcfile , status) IF_NOTOK_RETURN(status=1) call ReadRc(rcF, 'diffusion.write', diffusion_write, status, default=.false.) IF_ERROR_RETURN(status=1) if (diffusion_write) then call ReadRc(rcF, 'diffusion.dir', diffusion_dir, status) IF_NOTOK_RETURN(status=1) endif call ReadRc(rcF, 'my.levs', mlevs, status) IF_NOTOK_RETURN(status=1) call Done(rcF, status) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine Diffusion_Init ! *** subroutine Diffusion_Done( status ) ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Diffusion_Done' ! --- begin -------------------------------- ! nothing to be done ! ok status = 0 end subroutine Diffusion_Done ! ================================================================ !------------------------------------ ! ! Purpose: ! ------- ! this subroutine reads and prepares all fields needed for the calculation of kzz ! ! External ! -------- ! dd_get_3_hourly_surface_fields ! dd_calc_ustar_raero_rb ! dd_coarsen_fields ! ! Reference ! --------- ! Ganzeveld and Lelieveld (1996) and references therein. ! Adjusted for TM5, Bram Bregman, August 2003 ! ! subroutine calc_kzz( status ) use binas, only : grav use dims, only : im, jm, lmax_conv, idate, lm, okdebug use dims, only : nread, ndyn, tref, at, bt use dims, only : revert use global_data, only : mass_dat, wind_dat, region_dat #ifndef without_convection #ifndef without_diffusion use global_data, only : conv_dat #endif #endif use meteoData, only : spm_dat, pu_dat, pv_dat, m_dat use meteoData, only : temper_dat, humid_dat, gph_dat use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo, dist_arr_stat integer, intent(out) :: status ! --- local ------------------------------ character(len=10) :: c_time real,dimension(:,:,:),pointer :: phlb,gph,t,q,m,pu,pv real,dimension(:,:,:),pointer :: am,bm,cm ! fluxes in kg/timestep(region) #ifndef without_convection #ifndef without_diffusion real,dimension(:,:,:),pointer :: dkg !vert diff (ks/s) #endif #endif real,dimension(:),pointer :: dxyp integer,parameter :: avg_field=1, nlon360=360, nlat180=180 integer :: i,region, nsend, ntimes real, dimension(:,:,:),allocatable,target :: m_adj, phlb_adj real, dimension(:,:),allocatable :: p_adj integer :: imr, jmr, lmr, j,l real, allocatable, target :: phlb_t(:,:,:) real, allocatable, target :: m_t(:,:,:) integer :: i1, i2, j1, j2 character(len=*), parameter :: rname = mname//'/Calc_kzz' ! --- begin ---------------------------------- if (okdebug) write(*,*) 'start of calc_kzz' ! Allocate work arrays, and fill ghost cells DO region = 1, nregions CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ALLOCATE(ustar_loc(region)%surf( i1:i2, j1:j2)) ALLOCATE( sr_mix(region)%surf( i1:i2, j1:j2)) CALL UPDATE_HALO( dgrid(region), wind_dat(region)%am, wind_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), wind_dat(region)%bm, wind_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), wind_dat(region)%cm, wind_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ENDDO ! -- for output of time header on debug fields write(c_time,'(i4,3(i2))') idate(1:4) !CMK corrected do i=1,10 if (c_time(i:i)==' ') c_time(i:i)='0' enddo ! ! Surface data sets may have the following characteristics ! 1. instanteneous. The time written in the file is valid from t-1.5 to t+1.5 ! 2. Accumulated. The time written in the file is valid from t to t+3 ! 3. Daily averaged. The time on the file is always 00 hours, and valid until t+24 ! ! *** It is essential to understand this timing when reading and applying these sets. ! ! ! fd2mk it has to be decided to 'save' fields like vgrat_low and daily average surface fields ! for later use, or to read it every 3 hours time step REG: do region = 1,nregions ! local grid sizes lmr=lm(region) call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) if (revert == -1) then !mkadj !mass and pressure are at t=t+dt in adjoint !these should be brought back to t to get !exactly the same diffusion as in forward run allocate ( m_adj (i1:i2, j1:j2, lmr ) ) allocate ( phlb_adj(i1:i2, j1:j2, lmr+1) ) allocate ( p_adj (i1:i2, j1:j2 ) ) am => wind_dat(region)%am bm => wind_dat(region)%bm cm => wind_dat(region)%cm dxyp => region_dat(region)%dxyp ntimes = (nread/(ndyn))*tref(region) m_adj = m_dat(region)%data(i1:i2, j1:j2, 1:lmr) do i=1,ntimes m_adj = m_adj + revert*am(i1-1:i2-1, j1 :j2 , 1:lmr ) & !east - revert*am(i1 :i2 , j1 :j2 , 1:lmr ) & !west + revert*bm(i1 :i2 , j1 :j2 , 1:lmr ) & !south - revert*bm(i1 :i2 , j1+1:j2+1, 1:lmr ) & !north + revert*cm(i1 :i2 , j1 :j2 , 0:lmr-1) & !lower - revert*cm(i1 :i2 , j1 :j2 , 1:lmr ) !upper enddo p_adj(:,:) = 0.0 do l=1,lmr ! note that on the EDGES masses are coarse> ! diffusion iis applied only on core zoom> do j = j1, j2 do i = i1, i2 p_adj(i,j) = p_adj(i,j) + m_adj(i,j,l)*grav/dxyp(j) end do end do end do do l=1,lmr+1 do j = j1, j2 do i = i1, i2 phlb_adj(i,j,l) = at(l)+bt(l)*p_adj(i,j) end do end do end do deallocate(p_adj) nullify(am) nullify(bm) nullify(cm) nullify(dxyp) phlb => phlb_adj m => m_adj else ! Bug fix 2006-09-21 (AJS) ! The two half level pressure arrays phbl_t and phlb_k ! are sometimes not updated together when running in parallel ! with zoom regions. To avoid problems, pressure and mass ! are computed here given the surface pressure in the middle ! of the current time interval, stored in spm_dat. !-- !phlb => mass_dat(region)%phlb_t !m => mass_dat(region)%m_t !-- allocate( phlb_t(i1:i2, j1:j2, 1:lmr+1) ) do l = 1, lmr+1 phlb_t(:,:,l) = at(l) + bt(l)*spm_dat(region)%data(i1:i2, j1:j2, 1) end do allocate( m_t(i1-2:i2+2, j1-2:j2+2, 1:lmr) ) m_t = 0.0 do l = 1, lmr do j = j1, j2 m_t(i1:i2,j,l) = ( phlb_t(:,j,l) - phlb_t(:,j,l+1) )*region_dat(region)%dxyp(j)/grav end do end do phlb => phlb_t m => m_t !-- endif pu => pu_dat(region)%data pv => pv_dat(region)%data gph => gph_dat(region)%data t => temper_dat(region)%data q => humid_dat(region)%data #ifndef without_convection #ifndef without_diffusion dkg => conv_dat(region)%dkg #endif #endif ! ustar_loc call dd_calc_ustar(region, i1, i2, j1, j2) if (okdebug) then ! call dd_field_statistics(region, 'ustar_loc', ustar_loc(region)%surf, i1, j1, status) call dist_arr_stat(dgrid(region), 'ustar_loc', ustar_loc(region)%surf, 0, status) end if #ifndef without_convection #ifndef without_diffusion ! calculate kzz call bldiff( region, phlb, m, i1, i2, j1, j2 ) #endif #endif nullify(phlb) nullify(m) !-- deallocate( phlb_t ) deallocate( m_t ) !-- nullify(pu) nullify(pv) #ifndef without_convection #ifndef without_diffusion nullify(dkg) #endif #endif nullify(t) nullify(q) nullify(gph) if(revert == -1) then deallocate(m_adj) deallocate(phlb_adj) endif end do REG ! Done do region = 1,nregions deallocate(sr_mix(region)%surf) deallocate(ustar_loc(region)%surf) end do if (okdebug) write(*,*) 'end of calc_kzz' end subroutine calc_kzz !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: dd_calc_ustar ! ! !DESCRIPTION: ! ! Calculate friction velocity (ustar) ! ! Uses the log normal wind profile information stored by ECMWF in 10 meter wind ! to estimate a local ustar over land ! uses Charnock equation over sea to estimate ustar ! aerodynamic resistance is calculated from heat fluxes and ustar using a constant reference height ! sub laminar resistance from ustar ! ! Reference: ! o M.Z. Jacobson, "Fundamentals of Atmospheric Modeling", sections 8.3-8.4 ! ! !REVISION HISTORY: ! ! 2003 ??; Ge Verver (personal communication, 2003) ! ! 2010-08 Arjo Segers ! Debugged computation of u* over sea: ! abs. surface stress should be devided by air density. ! Set minum value for u* to trap division by zero (MK). ! ! 2012-01 P. Le Sager ! adapted for lon-lat MPI domain decomposition !------------------------ subroutine dd_calc_ustar(region, i1, i2, j1, j2) use binas, only : grav, vKarman use dims, only : okdebug use meteoData, only : lsmask_dat, sr_ecm_dat, sr_ols_dat use meteoData, only : wspd_dat use meteoData, only : slhf_dat, sshf_dat use meteoData, only : ewss_dat, nsss_dat !--- in/out ------------------------------------ integer,intent(in) :: region, i1, i2, j1, j2 ! --- const ------------------------------------- ! Garret relation: real,parameter :: alfa_garret = 0.11 ! air density at seal level and 293 K : real, parameter :: rho_air = 1.2 ! kg/m3 ! Charnock relation: real,parameter :: alfa_charnock = 0.018 ! kinematic viscosity of air at about 300 K : real,parameter :: nu_air = 1.5e-5 ! m2/s real, parameter :: Href=30. ! constant reference height for calculations of raero ! some constants specific to calculation of ustar and raero real,parameter :: rhoCp = 1231.0 real,parameter :: rhoLv = 3013000.0 real, parameter :: rz0 =2.0 !--- local ------------------------------------- integer :: i,j real :: tau_z_abs real :: ustar_sea,ustar_land,xland,sr_sea,sr_land,sr_help,wind10m ! --- begin ---------------------------------- do j=j1, j2 do i=i1, i2 xland=lsmask_dat(region)%data(i,j,1)/100. ! 0-1 ! SEA: ! ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling", ! section "8.3 Friction velocity" : ! The friction velocity : ! u* = sqrt( |tau_z|/rho_a ) [m/s] ! where: ! tau_z = (ewss,nsss) surface stress [N/m2=kg/m/s2] ! rho_air air density [kg/m3] ; about 1.2 at sea level and 293 K ! tau_z_abs = sqrt(ewss_dat(region)%data(i,j,1)**2+& nsss_dat(region)%data(i,j,1)**2 ) ! N/m2 ustar_sea = sqrt(tau_z_abs/rho_air) ! m/s ! ! limitation to avoid division by zero: ustar_sea = max( 0.001, ustar_sea ) ! m/s ; minium of 1 mm/s ! ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling", ! section "8.4 Surface roughness lengths" : ! "For smooth surfaces, such as over a smooth ocean with low wind speeds (Garret): ! z0m = 0.11 nu_a / u* ! where nu_a is the kinematic viscosity of air. ! Over a rough ocean with high wind speeds: ! z0m = alpha_c (u*)**2 / grav ! which is the 'Charnock relation', ! where alpha_c ~ 0.016 is the 'Charnock constant'." ! ! Here the sum of these two parameterizations is used, ! where the first part is dominant for u*<0.1 and the second for u*>0.1 ; ! minimum value for u* prevents division by zero: ! sr_sea = alfa_garret * nu_air / ustar_sea + & alfa_charnock * ustar_sea**2 / grav ! ! LAND ! calculate the 'local' surface roughness for momentum that is consistent with 10 m wind ! and the Olsson data base, we assume that the windspeed at 75 m is independent of ! surface roughness ! wind10m = wspd_dat(region)%data(i,j,1) if ( xland > 0.0 ) then !>>> TvN ! Over land, the friction velocity ustar ! is calculated from the 10 m wind speed, u10. ! In IFS u10 is a diagnostic variable, calculated to be ! compatible with "open-terrain" wind speed observations. ! Over rough or inhomogeneous terrain (z0M > 0.03 m), ! it is calculated using an aerodynamic roughness ! length typical for open terrain with grassland (0.03 m) ! (see IFS cy31r1 documentation, p. 46). ! In the expressions applied in TM5, ! the stability profile functions Psi_M have disappeared, ! and only the logarithmic terms are kept. ! Apart from this, the expression for ustar was incorrect: ! 10./sr_land should be 75./sr_land. ! This has now been corrected. ! Moreover, over islands and near coast lines, ! zeros can occur in sr_ols_dat, ! which have to be removed. ! However, a lower bound of 1e-2 m = 1 cm ! seems too high for smooth surfaces, ! like deserts and ice caps. ! The minimum positive value in sr_ols_dat ! is 2.5e-4 m = 0.025 cm, ! which is obtained in parts of Antarctica. ! This is likely the result of regridding of ! a field with minimum value of 0.1 cm ! from 0.5x0.5 to 1x1 degrees. ! Please note that with the introduction of CY31R1, ! the orographic contribution to the aerodynamic roughness length ! in IFS has been replaced by an explicit specification of stress ! on model levels due to turbulent orographic form drag, ! and the climatological variable previously used ! has been replaced by a prognostic variable. ! In TM5, the climatological variable is still used, ! but it would be better to use the prognostic variable instead. ! Note also that the same calculation ! is done in dry_deposition.F90. !sr_land=max(1e-2,sr_ols_dat(region)%data(i,j,1)) ! occurs at Islands, etc. sr_land=max(1e-3,sr_ols_dat(region)%data(i,j,1)) sr_help=min(sr_ecm_dat(region)%data(i,j,1),0.03) !fd ustar_land=vKarman*wind10m(i,j)/alog(10./sr_help) !ustar consistent with 'clipped' large scale roughness !ustar_land=vKarman*wind10m/& ! alog(10./sr_help)*alog(75./sr_help)/alog(10./sr_land) ustar_land=vKarman*wind10m/& alog(10./sr_help)*alog(75./sr_help)/alog(75./sr_land) ! <<< TvN else sr_land = 0.0 ustar_land = 0.0 endif ! interpolate between land and sea for mixed pixels: ustar_loc(region)%surf(i,j) = xland*ustar_land+(1-xland)*ustar_sea sr_mix (region)%surf(i,j) = xland* sr_land+(1-xland)* sr_sea end do !i end do !j if (okdebug) write(*,*) 'end calc_ustar' end subroutine dd_calc_ustar ! *** #ifndef without_convection #ifndef without_diffusion !------------------------------------------------------------------ ! Procedure: ! (1) calc_Kv: Diffusion coefficients full atmosphere from shear ! (2) pblhght: Determine the boundary layer height ! (3) difcoef: Replace diffusion coefficients in the boundary layer !------------------------------------------------------------------ ! subroutine adjusted for TM5, Bram Bregman, August 2003. !-------------------------IO--------------------------------------- subroutine bldiff( region, phlb, m, i1, i2, j1, j2 ) use binas, only: ae, cp_air, Rgas, grav, Lvap use binas, only: xmair use binas, only: vkarman use dims, only: im, jm, lm, lmax_conv use dims, only: okdebug use dims, only: gtor, dx, dy, ybeg, xref, yref use dims, only: isr, ier, jsr, jer use toolbox, only: dumpfield use global_data, only : conv_dat use meteodata , only : slhf_dat, sshf_dat, pu_dat, pv_dat, gph_dat, temper_dat, humid_dat ! --- in/out ---------------------------------------------- integer,intent(in) :: region, i1, i2, j1, j2 real,dimension(:,:,:),pointer :: phlb, gph, t, q, m, pu, pv, dkg ! --- const ---------------------------------- ! minimum value diffusion coefficient real, parameter :: ckmin = 1.e-15 real, parameter :: sffrac= 0.1 real, parameter :: onet = 1./3. real, parameter :: rair = Rgas*1.e3/xmair real, parameter :: ccpq = 1869.46/cp_air -1. real, parameter :: epsilo = rair/461.51 real, parameter :: apzero = 100000. ! reference pressure (usually 1000 hPa) real, parameter :: ccon = sffrac*vkarman real, parameter :: betam=15.,betah=15. real, parameter :: binm = betam*sffrac real, parameter :: binh = betah*sffrac real, parameter :: fr_excess = 1. real, parameter :: fak = fr_excess * 8.5 real, parameter :: fakn = fr_excess * 7.2 real, parameter :: zkappa = rair / cp_air real, parameter :: zcrdq = 1.0/epsilo - 1.0 real, parameter :: tiny = 1.e-9 ! bound wind**2 to avoid dividing by 0 real, parameter :: zacb = 100. !factor in ustr-shearproduction term real, parameter :: ricr = 0.3 ! critical richardson number real, parameter :: ssfrac=0.1 ! --- local ------------------------------------ real, dimension(:,:,:), allocatable:: wkvf,& ! outer layer Kv (at the top of each layer) zdup2,& ! vertical shear squared zrinub,& ! Richardson number wthm,& ! potential temperature zthvk,& ! virtual potential temperature pf,& ! full level pressure whm, & ! height of layer centers uwind,vwind,& ! wind velocities [m s-1] kvh ! real, dimension(:,:), allocatable :: ustr,& ! velocity scale wheat,& ! surface sensible heat flux [K m/s] wqflx,& ! surface latent heat flux [kg m/ (kg s)] wobkl,& ! Monin Obukhov length wheatv,& ! surface virtual heat flux [K m/s] wtseff ! theta_v at lowest level [K] integer :: i, j, l real :: intfac !cmk BUG was integer! ! variables for the calculation of free atmosphere vertical diffusion coefficients wkvf real :: cml2,zlambdac,arg,zdz,zthva,zfunst, zfstabh,zsstab,zkvn,dyy real, dimension(j1:j2) :: dxx, lat integer :: jq! jqif function ! variabales for the calculation of the boundary layer height real :: wrino ! Richardson number real :: thvparc ! parcel virtual potential temperature real :: fmt,wsc,tkv,zrino,vvk,dtkv,zexcess integer :: jwork real,dimension(:,:),pointer :: pblh ! planetary boundary layer height [m] ! variables to calculate difusion coefficient in the boundary layer !real :: wsc,cml2 real :: z,zwstr,zkvh,kvhmin,zm,zp,zslask,zsl1, & zh,zzh,zl,pr,zpblk,zfstabh1,zfstabh2,& term1,term2,term3,obukov integer :: jq1,jq2,jq3,jq4,jq5,jck,jqq, jq6, jq7 real :: yres real :: xres real :: thvgrad, kvhentr ! --- begin ------------------------------------------ if (okdebug) write(*,*) 'start of bldiff' allocate(wkvf ( i1:i2, j1:j2, lm(region))) allocate(zdup2 ( i1:i2, j1:j2, lm(region))) allocate(zrinub( i1:i2, j1:j2, lm(region))) allocate(wthm ( i1:i2, j1:j2, lm(region))) allocate(zthvk ( i1:i2, j1:j2, lm(region))) allocate(pf ( i1:i2, j1:j2, lm(region))) allocate(whm ( i1:i2, j1:j2, lm(region))) allocate(uwind ( i1:i2, j1:j2, lm(region))) allocate(vwind ( i1:i2, j1:j2, lm(region))) allocate(kvh ( i1:i2, j1:j2, lm(region))) allocate(ustr ( i1:i2, j1:j2 )) allocate(wheat ( i1:i2, j1:j2 )) allocate(wqflx ( i1:i2, j1:j2 )) allocate(wobkl ( i1:i2, j1:j2 )) allocate(wheatv( i1:i2, j1:j2 )) allocate(wtseff( i1:i2, j1:j2 )) pblh => conv_dat(region)%blh pu => pu_dat(region)%data pv => pv_dat(region)%data gph => gph_dat(region)%data t => temper_dat(region)%data q => humid_dat(region)%data dkg => conv_dat(region)%dkg ! mid-layer pressure levels ! potential temperature field ! virtual potential temperature field ! !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (lm, jsr, jer, isr, ier, region, pf, phlb, t, q, wthm, & !PLS $OMP zthvk, whm, gph) & !PLS $OMP private (i, j, l, jss, jee, intfac) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region) do j=j1,j2 do i=i1,i2 pf(i,j,l) = (phlb(i,j,l) + phlb(i,j,l+1) )*0.5 wthm(i,j,l) = t(i,j,l) * ( apzero/pf(i,j,l) )**zkappa zthvk(i,j,l) = wthm(i,j,l) * (1.0 + zcrdq*q(i,j,l)) ! mid-layer heights (note: gph array is for levels 1:lm+1) intfac = (phlb(i,j,l)-pf(i,j,l)) / (phlb(i,j,l)-phlb(i,j,l+1)) ! cmk: intfac is now real! whm(i,j,l) = (gph(i,j,l)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1) -gph(i,j,1)) * intfac ! cmk bug corrected: ! cmk old: whm(i,j,l) = (gph(i,j,l+1)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1) -gph(i,j,1)) * intfac end do end do end do !PLS $OMP END PARALLEL ! convert units of sensible and latent heat fluxes ! sshf: W/m2 / (J/kg/K) / (kg/m3) = J/s/m2 /J kg K /kg m3 = K m/s ! slhf: W/m2 / (J/kg) / (kg/m3) = J/s/m2 /J kg /kg m3 = m/s ! ustr is local velocity scale, NOT from ECMWF surface stresses ! !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (jm, im, region, wheat, sshf_dat, gph, phlb, slhf_dat, & !PLS $OMP wqflx, ustr, ustar_loc) & !PLS $OMP private (i, j, js, je) ! call my_omp_domain (1, jm(region), js, je) do j = j1,j2 do i = i1,i2 wheat(i,j) = -sshf_dat(region)%data(i,j,1) * (gph(i,j,2) - & gph(i,j,1)) * grav / (phlb(i,j,1) - & phlb(i,j,2)) / cp_air wqflx(i,j) = -slhf_dat(region)%data(i,j,1) * (gph(i,j,2) - & gph(i,j,1)) * grav / (phlb(i,j,1) - & phlb(i,j,2)) / lvap ustr(i,j) = max (ustar_loc(region)%surf(i,j), 0.01) end do end do !PLS $OMP END PARALLEL ! !----------------------------------------------------------------- ! compute the free atmosphere vertical diffusion coefficients wkvf ! height-dependent asymptotic mixing length implemented ! Holtslag and Boville, 1993, J. Climate, 6, 1825-1842. !------------------------------------------------------------------ ! first calculate wind velocities [m s-1] from the mass fluxes yres = dy/yref(region) xres = dx/xref(region) do j = j1,j2 lat(j) = ybeg(region) + 0.5 * yres + yres * (j-1) enddo dxx(j1:j2) = ae * xres * gtor * cos(lat(j1:j2)*gtor) dyy = ae * yres * gtor !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (lm, jsr, jer, isr, ier, region, dxx, dyy, pu, pv, m, & !PLS $OMP uwind, vwind) & !PLS $OMP private (i, j, l, jss, jee) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region) do j=j1,j2 do i=i1,i2 uwind(i,j,l) = dxx(j)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l) vwind(i,j,l) = dyy *(pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l) enddo enddo enddo !PLS $OMP END PARALLEL ! ! Initialize gradient Richardson number and free tropospheric k values ! if (okdebug) write(*,*) ' Richardson number' !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (lm, jsr, jer, isr, ier, region, gph, zdup2, uwind, & !PLS $OMP vwind, whm, phlb, pf, zthvk, zrinub, wkvf) & !PLS $OMP private (i, j, l, jss, jee, jq, zlambdac, cml2, zdz, arg, & !PLS $OMP zthva, zsstab, zfunst, zfstabh, zkvn) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region)-1 do j=j1,j2 do i=i1,i2 ! ! if (delta gph < 1000) then jq=1 jq = jqif( 1000.0, (gph(i,j,l+1)-gph(i,j,1)) ) zlambdac=jq*300.+(1-jq)*(30. + 270. * exp (1.-(gph(i,j,l+1)-gph(i,j,1))/1000.)) cml2=(1./( 1./zlambdac + 1./(vkarman*(gph(i,j,l+1)-gph(i,j,1))) ))**2. ! ! vertical shear squared, ! min value of (delta v)**2 prevents zero shear ! zdup2(i,j,l) = (uwind(i,j,l+1)-uwind(i,j,l))**2 + & (vwind(i,j,l+1)-vwind(i,j,l))**2 zdup2(i,j,l) = max(zdup2(i,j,l),1.e-10) zdz = whm(i,j,l+1) - whm(i,j,l) zdup2(i,j,l) = zdup2(i,j,l)/(zdz**2) ! ! static stability (use virtual potential temperature) ! dv now calculated at interface height (lineair in pressure) ! arg=(log(phlb(i,j,l+1))-log(pf(i,j,l)))/(log(pf(i,j,l+1))-log(pf(i,j,l))) zthva = zthvk(i,j,l) + arg * (zthvk(i,j,l+1)-zthvk(i,j,l)) zsstab = grav/zthva*(zthvk(i,j,l+1) - zthvk(i,j,l))/zdz ! ! gradient Richardson number ! zrinub(i,j,l) = zsstab/zdup2(i,j,l) ! ! stability functions ! zfunst = max(1. - 18.*zrinub(i,j,l),0.) zfunst = sqrt(zfunst) ! zfstabh = 1.0 / (1.0 + 10.0*zrinub(i,j,l)*(1.0+8.0*zrinub(i,j,l))) ! ! neutral diffusion coefficient ! zkvn = cml2 * sqrt(zdup2(i,j,l)) ! ! correction with stability functions ! jq = jqif( zrinub(i,j,l),0.0 ) ! thus: if (zrinub>0) then jq=1 wkvf(i,j,l)=jq*(max(ckmin,zkvn*zfstabh))+(1-jq)*(max(ckmin,zkvn*zfunst)) end do end do end do !PLS $OMP END PARALLEL ! ! Determine the boundary layer height !-------------------------------------------------------------------------------- ! based on: Vogelezang and Holtslag, 1996, Bound. Layer Meteorol., 81, 245-269. !-------------------------------------------------------------------------------- if (okdebug) write(*,*) ' boundary layer height' !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (lm, jsr, jer, isr, ier, region, wtseff, wthm, q, & !PLS $OMP wheatv, wheat, wqflx, wobkl, ustr, pblh, uwind, vwind, & !PLS $OMP whm ) & !PLS $OMP private (i, j, l, jss, jee, jwork, wrino, vvk, tkv, dtkv, & !PLS $OMP zrino, jq, fmt, wsc, zexcess, thvparc ) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do j=j1,j2 do i=i1,i2 ! compute bottom level virtual temperature and heat flux and Monin-Obukhov Length wtseff(i,j) = wthm(i,j,1)*(1. + zcrdq*q(i,j,1)) wheatv(i,j) = wheat(i,j) + zcrdq*wthm(i,j,1)*wqflx(i,j) wobkl(i,j) = -wtseff(i,j)*ustr(i,j)**3 / (grav*vkarman*(wheatv(i,j)+sign(1.e-10,wheatv(i,j))) ) ! initialize pblh(i,j) = 0.0 jwork = 1 wrino = 0.0 do l=1,lm(region) vvk = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + & zacb*ustr(i,j)*ustr(i,j) vvk = max(vvk,tiny) tkv = wthm(i,j,l)*(1. + zcrdq*q(i,j,l)) dtkv = tkv-wtseff(i,j) ! ! Bulk Richardson number ! zrino = grav*dtkv * (whm(i,j,l)-whm(i,j,1)) / (wtseff(i,j)*vvk) zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero zrino = jwork*zrino jq = jqif(ricr,zrino) ! thus: if (zrino < ricr) then jq=1 ! ! calculate pblh from linear interpolation in Ri(bulk) ! if ( l == 1 ) then pblh(i,j) = jq*pblh(i,j) + (1-jq) * & ( (ricr-wrino)/(zrino-wrino)*(whm(i,j,l) ) ) else pblh(i,j) = jq*pblh(i,j) + (1-jq) * & ( whm(i,j,l-1) + (ricr-wrino)/(zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1)) ) end if ! ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0 ! and pblh does not change anymore ! jwork = jq*jwork ! ! if pblh is found already avoid dividing by zero next level by a faked value ! of wrino (jwork=0 already) ! wrino = zrino+(1-jwork)*0.1 end do pblh(i,j) = (1-jwork)*pblh(i,j) + whm(i,j,lm(region))*jwork ! ! second calculation of pblh including excess surface temperature using ! a convective velocity scale wsc (wm, not equal to w*!!!), using result ! of the first calculation of pblh ! jq = jqif(wheatv(i,j),0.) ! thus: if (wheatv(i,j) > 0.) then jq=1 jwork = jq fmt = ( jq*(1.0 - binm*pblh(i,j)/wobkl(i,j)) + (1.0-jq) )**onet wsc = ustr(i,j)*fmt zexcess = wheatv(i,j)*fak/wsc ! ! improve pblh estimate under convective conditions using ! convective temperature excess (zexcess) ! thvparc = wtseff(i,j)+jq*zexcess ! jwork = 1 wrino = 0.0 ! do l=2,lm(region) ! ! Bulk Richardson number: ! if stable corrected with extra shear term ! if unstable now also corrected for temperature excess ! vvk = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + & zacb*ustr(i,j)*ustr(i,j) vvk = max(vvk,tiny) tkv = wthm(i,j,l)*(1. + zcrdq*q(i,j,l)) zrino = grav*(tkv-thvparc) * (whm(i,j,l)-whm(i,j,1))/(vvk*wtseff(i,j)) zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero zrino = zrino*jwork jq = jqif(ricr,zrino) ! thus: if (zrino < ricr) then jq=0 ! ! calculate pblh from linear interpolation in Ri(bulk) ! pblh(i,j) = jq*pblh(i,j)+(1-jq)* (whm(i,j,l-1)+(ricr-wrino)/ & (zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1))) ! ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0 ! and pblh does not change anymore ! jwork = jq*jwork ! ! if pblh is found already avoid dividing by zero next level by a faked value ! of wrino (jwork=0 already) ! wrino = zrino+(1-jwork)*0.1 end do pblh(i,j) = (1-jwork)*pblh(i,j) + jwork*whm(i,j,lm(region)) pblh(i,j) = max(pblh(i,j), 100.0) ! set minimum value for pblh end do end do !PLS $OMP END PARALLEL ! ! Diffusion coefficients in the surface layer ! !-------------------------------------------------------------------------------- ! the revised LTG scheme (Beljaars and Viterbo, 1998) ! Modification - April 1, 1999: prevent zpblk to become too small in very stable conditions !-------------------------------------------------------------------------------- if (okdebug) write(*,*) ' diffusion coeff' ! initialize output array with minimum value kvhmin = 0.1 ! ! Loop to calculate the diffusivity ! !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (lm, jsr, jer, isr, ier, region, pblh, gph, wobkl, & !PLS $OMP wheatv, whm, zrinub, zdup2, kvhmin, wkvf, ustr, & !PLS $OMP wtseff, kvh, zthvk) & !PLS $OMP private (i, j, l, jss, jee, jq, z, zh, zl, zzh, jq1, jq2, jq3, & !PLS $OMP jq4, jq5, jq6, jq7, jqq, cml2, zfstabh, zpblk, zkvh, & !PLS $OMP zslask, fmt, wsc, zwstr, obukov, term1, term2, term3, & !PLS $OMP pr, thvgrad, kvhentr) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region)-1 do j=j1,j2 do i=i1,i2 jq = jqif(pblh(i,j),(gph(i,j,l+1)-gph(i,j,1))) ! if top of level hh < pblh jq=1 inside BL z = gph(i,j,l+1)-gph(i,j,1) zh = jq*z/pblh(i,j) ! only defined in BL (jq=1) ! zl must not be zero zl = jq*z/wobkl(i,j)+(1-jq) ! z/L outside BL = 1 zzh = jq*(1.-zh)**2 ! (1-z/h)^2 only in BL jq1 = jqif(wheatv(i,j),0.) ! jq1 is 1 in unstable case jq2 = jq*(1-jq1) ! jq2 is 1 in stable BL jq3 = jq*jq1 ! jq3 is 1 in unstable BL jq6 = jqif(pblh(i,j), whm(i,j,l)) jq7 = jqif(whm(i,j,l+1), pblh(i,j)) jq6 = jq1*jq6*jq7 ! jq6 is 1 at at the kvh-level that is located between the ! mid-levels that surround pblh jq1 = jqif(sffrac,zh) ! jq1 is 1 in surface layer jq4 = jq3*jq1 ! jq4 is 1 in unstable surface layer jq5 = jq3*(1-jq1) ! jq5 is 1 outside unstable surface layer ! ! calculate coefficients for momentum, the values for heat are obtained by dividing by Prantl number ! ! stable and neutral: ! cml2 = (1./( 1./(vkarman*z) + 1./450. ))**2. jqq = jqif(zrinub(i,j,l),0.) ! ! stability function for momentum ! zfstabh = jqq/(1.+10.*zrinub(i,j,l)*sqrt(1.+jqq*zrinub(i,j,l))) + (1-jqq) zpblk = cml2 * sqrt(zdup2(i,j,l))*zfstabh zkvh = jq2*max(zpblk,kvhmin)+(1-jq2)*wkvf(i,j,l) ! take free troposp. values outside BL ! ! unstable case in surface layer ! zslask = 1.-betah*zl zslask = jq4*zslask+(1-jq4) zpblk = (1-jq4)*zkvh+jq4*(ustr(i,j)*vkarman*z*zzh*zslask**(1./3.)) ! ! unstable case above surface layer ! ! evaluate stability function for momentum at height z = 0.1*pblh (top of surface layer) ! zslask = 1.-ssfrac*betah*pblh(i,j)/wobkl(i,j) zslask = jq5*zslask+(1-jq5) fmt = zslask**(1./3.) ! ! use surface layer definition for w_m wsc = ustr(i,j)*fmt zpblk = (1-jq5)*zpblk + jq5*wsc*vkarman*z*zzh ! ! Determine Prandt Number ! ! Pr-number is assumed to be constant throughout the CBL, i.e. in surface and mixed layer ! NOTE: this is different from the original formulation ! zwstr not used if wheatv<0. ! zwstr = (abs(wheatv(i,j))*grav*pblh(i,j)/wtseff(i,j))**(1./3.) ! bh obukov = jq3*wobkl(i,j)-(1-jq3) term1 = (1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.) term2 = sqrt(1.-ssfrac*betah*pblh(i,j)/obukov) term3 = ccon*fakn*zwstr/(ustr(i,j)*(1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.)) pr = jq3*(term1/term2+term3)+(1-jq3) ! ! NOTE that kvh applies to the top of the level ! kvh(i,j,l) = jq3*max(zpblk/pr,kvhmin)+(1-jq3)*zkvh ! if in the entrainment zone, override with prescribed entrainment thvgrad = (zthvk(i,j,l+1)-zthvk(i,j,l))/(whm(i,j,l+1)-whm(i,j,l)) kvhentr = 0.2*wheatv(i,j)/thvgrad kvh(i,j,l) = jq6*kvhentr + (1-jq6)*kvh(i,j,l) end do end do end do !PLS $OMP END PARALLEl ! calculate vertical diffusion !PLS $OMP PARALLEL & !PLS $OMP default (none) & !PLS $OMP shared (jsr, jer, isr, ier, region, dkg, kvh, m, gph) & !PLS $OMP private (i, j, l, jss, jee) ! call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lmax_conv-1 do j=j1,j2 do i=i1, i2 dkg(i,j,l)=max(0.,kvh(i,j,l))*2.*(m(i,j,l+1)+m(i,j,l))/(gph(i,j,l+2)-gph(i,j,l))**2 ! CMK dec 2002 ----> gph changed to 1--->lm+1 boundaries enddo enddo enddo dkg(i1:i2,j1:j2,lmax_conv) = 0.0 !PLS $OMP END PARALLEL deallocate(wkvf) deallocate(zdup2) deallocate(zrinub) deallocate(wthm) deallocate(zthvk) deallocate(pf) deallocate(whm) deallocate(ustr) deallocate(wheat) deallocate(wqflx) deallocate(wobkl) deallocate(wheatv) deallocate(wtseff) deallocate(uwind) deallocate(vwind) deallocate(kvh) nullify(pblh) if (okdebug) write(*,*) 'end of bldiff' end subroutine bldiff #endif /* diffusion */ #endif /* convection */ ! ! define vectorizable 'if' function ! if xxxz>yyyz jqif=1 else jqif=0 ! INTEGER FUNCTION JQIF( xxxz, yyyz ) ! --- in/out --------------------------- real, intent(in) :: xxxz, yyyz ! --- begin ---------------------------- ! !jqif = nint( 0.5 - sign( 0.5, -dim(xxxz,yyyz)) ) ! jqif = nint( 0.5 + sign( 0.5, xxxz-yyyz ) ) end function jqif !===================================================================================================== !===================================================================================================== subroutine diffusion_filename(tau, region, fname) use datetime, only : tau2date use Go, only : pathsep use dims, only : region_name, revert integer(kind=8), intent(in) :: tau integer, intent(in) :: region character(len=1000), intent(out) :: fname character(len=*), parameter :: rname = mname//'/diffusion_filename' integer, dimension(6) :: idater if ( revert == 1 ) then call tau2date( tau, idater ) else ! read field from 3 hours before call tau2date( tau - 3600*3, idater ) end if write(fname, '(6a, i4.4, a1, i2.2, a1, i2.2, a1, "dkg_", i4, 4i2.2, ".nc")') & trim(diffusion_dir), pathsep, trim(mlevs), pathsep, region_name(region), pathsep, & idater(1), pathsep, idater(2), pathsep, idater(3), pathsep, idater(1:5) end subroutine diffusion_filename !===================================================================================================== !===================================================================================================== subroutine read_diffusion( status ) ! ! Read vertical diffusion from file and store in conv_dat%dkg ! ! Output: ! o status 0 = ok ! -1 = file not available for at least one region ! 1 = read error ! ! 19 Dec 2012 - P. Le Sager - TM5-MP version, using MDF ! --- modules ------------------------------ use dims, only : itau, revert, region_name, lm, okdebug use global_data, only : conv_dat use datetime, only : tau2date use MDF, only : MDF_OPEN, MDF_NETCDF, MDF_Inq_VarID, MDF_Get_Var, MDF_Close, MDF_READ use Partools, only : isRoot, par_broadcast use meteodata, only : global_lli use tm5_distgrid, only : dgrid, scatter ! --- in/out ---------------------------------------------- integer, intent(out) :: status ! --- const ----------------------------------------------- character(len=*), parameter :: rname = mname//', read_diffusion' ! --- local ----------------------------------------------- character(len=1000) :: fname integer :: region, fid, varid, imr, jmr, sz(3) integer, dimension(6) :: idater logical :: exist real, allocatable :: global_arr(:,:,:) ! --- begin ----------------------------------------------- if ( revert == 1 ) then call tau2date( itau, idater ) else ! read field from 3 hours before call tau2date( itau - 3600*3, idater ) end if do region = 1, nregions ! entire region grid size imr = global_lli(region)%nlon jmr = global_lli(region)%nlat ! for gathering sz = shape(conv_dat(region)%dkg) ! to get 3rd dimension if(isRoot)then allocate(global_arr(imr,jmr,sz(3))) else allocate(global_arr(1,1,1)) end if ! Create file name call diffusion_filename(itau,region,fname) if (isRoot) then inquire( file=fname, exist=exist ) end if call par_broadcast(exist, status) if ( .not. exist ) then write(gol,'(a,": dkg file """,a,""" not found, will create.")') rname, trim(fname) call goPr status = -1 return end if if (isRoot) then ! Open file ! This is a bit too verbose, leading to long output files. ! write(*,'(a,": reading ", a)') rname, trim(fname) call MDF_OPEN(TRIM(fname), MDF_NETCDF, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) ! Read dkg CALL MDF_Inq_VarID( fid, 'dkg', varid, status ) IF_ERROR_RETURN(status=1) CALL MDF_Get_Var( fid, varid, global_arr, status ) IF_NOTOK_RETURN(status=1) end if call scatter( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status) IF_NOTOK_RETURN(status=1) if(isRoot) then ! Read blh CALL MDF_Inq_VarID( fid, 'blh', varid, status ) IF_ERROR_RETURN(status=1) CALL MDF_Get_Var( fid, varid, global_arr(:,:,1), status ) IF_NOTOK_RETURN(status=1) end if call scatter( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) ! Close file if (isRoot) then call MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) end if deallocate(global_arr) end do status = 0 end subroutine read_diffusion !===================================================================================================== !===================================================================================================== subroutine write_diffusion( status ) ! ! Write vertical diffusion to files (for all regions) ! ! Output: ! o status 0 = ok, else not ok ! ! 19 Dec 2012 - P. Le Sager - made a TM5-MP version , that works using MDF ! --- modules ------------------------------ use dims, only : itau, revert, region_name, lm, okdebug use global_data, only : conv_dat use datetime, only : tau2date,date2tau use MDF, only : MDF_CREATE, MDF_NETCDF4, MDF_Def_Dim, MDF_Def_Var, MDF_Close use MDF, only : MDF_DOUBLE, MDF_NEW, MDF_Put_Var, MDF_Put_Att, MDF_EndDef, MDF_UNLIMITED use Partools, only : isRoot use meteodata, only : global_lli use tm5_distgrid, only : dgrid, gather use GO, only : pathsep ! --- in/out ---------------------------------------------- integer, intent(out) :: status ! --- const ----------------------------------------------- character(len=*), parameter :: rname = mname//'write_diffusion' ! --- local ----------------------------------------------- character(len=1000) :: fname, targetdir integer :: region, fid, dkgid, blhid, dimlon, dimlat, dimtime, dimlev, imr, jmr, sz(3) integer :: lonid, latid, timeid integer, dimension(6) :: idater real, allocatable :: global_arr(:,:,:) real :: sec_1970 integer(kind=8) :: itau_1970 integer :: ipos ! --- begin ----------------------------------------------- if ( revert == 1 ) then call tau2date( itau, idater ) else ! read field from 3 hours before call tau2date( itau - 3600*3, idater ) end if status=0 do region = 1, nregions ! entire region grid size imr = global_lli(region)%nlon jmr = global_lli(region)%nlat ! for gathering sz = shape(conv_dat(region)%dkg) ! to get 3rd dimension if(isRoot)then allocate(global_arr(imr,jmr,sz(3))) else allocate(global_arr(1,1,1)) end if ! Create and define if(isRoot) then ! Create file name call diffusion_filename(itau,region,fname) ! get the target directory from fname ipos = scan(trim(fname), pathsep, back=.true.) targetdir = fname(1:ipos-1) call system('mkdir -p ' // targetdir) if(okdebug) then write(gol,'(a,": creating ", a)') rname, trim(fname); call goPr end if call MDF_CREATE(TRIM(fname), MDF_NETCDF4, MDF_NEW, fid, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Dim( fid, 'longitude', imr, dimlon, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Dim( fid, 'latitude', jmr, dimlat, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Dim( fid, 'time', MDF_UNLIMITED, dimtime, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Var( fid, 'latitude', MDF_DOUBLE, (/dimlat/), latid, status=status) IF_NOTOK_RETURN(status=1) call MDF_Def_Var( fid, 'longitude', MDF_DOUBLE, (/dimlon/), lonid, status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, lonid, "units", values = "degrees_east", status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, latid, "units", values = "degrees_north", status=status) IF_NOTOK_RETURN(status=1) call MDF_Def_Var(fid, 'time', MDF_DOUBLE, (/dimtime/), timeid, status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, timeid, "units", values="seconds since 1970-01-01 00:00:00", status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, timeid, "comment", values="time values represent the beginning of the period", status=status) IF_NOTOK_RETURN(status=1) call MDF_Def_Dim( fid, 'levels', sz(3), dimlev, status ) IF_NOTOK_RETURN(status=1) call MDF_Def_Var( fid, 'dkg', MDF_DOUBLE, (/dimlon,dimlat,dimlev,dimtime/), dkgid, & status=status, deflate_level=3) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, dkgid, "units", values = "ks s-1", status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, dkgid, "long_name", values="vertical diffusivity coefficient", status=status) IF_NOTOK_RETURN(status=1) call MDF_Def_Var( fid, 'blh', MDF_DOUBLE, (/dimlon,dimlat,dimtime/), blhid, & status=status, deflate_level=3) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, blhid, "units", values = "m", status=status) IF_NOTOK_RETURN(status=1) call MDF_Put_Att(fid, blhid, "long_name", values="boundary layer height", status=status) IF_NOTOK_RETURN(status=1) call MDF_EndDef(fid, status) end if ! gather and write dkg/blh call gather( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot) then call MDF_Put_Var( fid, latid, global_lli(region)%lat_deg, status) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( fid, lonid, global_lli(region)%lon_deg, status) IF_NOTOK_RETURN(status=1) call date2tau([1970, 1, 1, 0, 0, 0], itau_1970) sec_1970 = real(itau) - real(itau_1970) call MDF_Put_Var(fid, timeid, (/sec_1970/), status) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( fid, dkgid, global_arr, status) IF_NOTOK_RETURN(status=1) end if call gather( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) if(isRoot) then call MDF_Put_Var( fid, blhid, global_arr(:,:,1), status) IF_NOTOK_RETURN(status=1) ! Close file call MDF_CLOSE( fid, status ) IF_NOTOK_RETURN(status=1) end if deallocate(global_arr) end do status = 0 end subroutine write_diffusion !===================================================================================================== !===================================================================================================== END MODULE DIFFUSION