!############################################################################### ! !ProTeX: 1.14-AJS ! !BOI ! ! !TITLE: TMM - TM Meteo ! !AUTHORS: Arjo Segers ! !AFFILIATION: KNMI ! !DATE: 21/04/2004 ! ! !INTRODUCTION: Introduction ! ! Module to access TM meteo. ! ! The main structure provides access to a list of opened meteo files. ! If a new meteo field is required, the subroutines search in the ! list wether the field is available. ! If not, a new file is opened and added to the list. ! Optionally, a shell script is invoked to search for a file ! and to store it locally if necessary. ! ! ! !INTRODUCTION: Usage ! ! \bv ! ! ! --- modules ----------------------------------------- ! ! use GO, only : TDate, NewDate ! ! use grid, only : TllGridInfo, Init, Done ! use grid, only : TLevelInfo, Init, Done ! ! use TMM, only : TTmMeteo ! use TMM, only : Init, Done ! use TMM, only : ReadField, ReadUVSP ! ! ! --- local ----------------------------------------- ! ! type(TllGridInfo) :: lli ! type(TLevelInfo) :: levi_ec, levi ! ! type(TTmMeteo) :: tmmd ! ! type(TDate) :: tday, t1, t2 ! ! real :: psurf(120,90) ! real :: temper(120,90,25) ! real :: pu(0:120,90,25) ! ! ! --- begin ----------------------------------------- ! ! ! define horizontal grid ! call Init( lli, -178.5, 3.0, 120, -89.0, 2.0, 90, status ) ! ! ! define vertical hybride levels ! call Init( levi_ec, 'ec60', status ) ! ecmwf levels ! call Init( levi, levi_ec, (/60,..,0/), status ) ! tm half level selection ! ! ! setup TM meteo access: ! call Init( tmmd, 'tm5.rc', status ) ! ! ! define time range for field: ! tday = NewDate( year=2003, month=01, day=01 ) ! t1 = NewDate( year=2002, month=12, day=31, hour=21 ) ! t2 = NewDate( year=2003, month=01, day=01, hour=03 ) ! ! ! Read meteo arrays; specify grid, parameter, time, etc. ! ! ! ! Type of grid is defined by nuv key: ! ! 'n' = normal grid (cell centers) ! ! 'u' = u-grid (east/west boundaries) ! ! 'v' = v-grid (north/south boundaries) ! ! ! ! Requested grid (lli/levi) might be different from the grid in the file. ! ! Horizontal: ! ! o if file contains same resolutions as defined by lli, ! ! a part of the data in the file is selected; ! ! o if file contains higher resolution fields, ! ! the ll array is filled with combined values from the file; ! ! depending on the parameter, the result is summed/avaraged/etc. ! ! Vertical: ! ! o if a file contains a supperset of the levels in levi, ! ! some levels are combined; ! ! depending on the parameter, the result is summed/avaraged/etc; ! ! o the file might contain fields with reversed level order. ! ! ! ! 3D fields require surface pressure for the level definition. ! ! It should be valid for [t1,t2] ! ! ! Best is to use the spm from ReadUVSP. ! ! call ReadUVSP ( tmmd, 'tmpp:od-fc-ml60-glb3x2', tday, t1, t2, lli, levi, sp1, spm, sp2, pu, pv, status ) ! ! call ReadField( tmmd, 'tmpp:od-fc-ml60-glb3x2', 'T' , tday, t1, t2, lli, 'n', levi, spm, temper, status ) ! call ReadField( tmmd, 'tmpp:od-fc-ml60-glb3x2', 'pu', tday, t1, t2, lli, 'u', levi, spm, pu , status ) ! ! ! ! ! TMPP surface fields can be called using the 1x1 as well as the 3x2 key. ! ! ! call ReadField( tmmd, 'tmpp:od-fc-sfc-glb1x1' , 'ci', tday, t1, t2, lli, 'n', ci, status ) ! call ReadField( tmmd, 'tmpp:od-fc-ml60-glb3x2', 'ci', tday, t1, t2, lli, 'n', ci, status ) ! ! ! ! Similar for spm surface pressures: ! ! ! call ReadField( tmmd, 'tmpp:od-fc-ml1-glb3x2' , 'spm', tday, t1, t1, lli, 'n', spm, status ) ! call ReadField( tmmd, 'tmpp:od-fc-ml60-glb3x2', 'spm', tday, t1, t1, lli, 'n', spm, status ) ! ! ! ! *** output ! ! call WriteField( tmmd, 'tmpp:od-fc-ml60-glb3x2', 'T', 'K', t1, t2, & ! lli, 'n', levi, spm, temper, status ) ! ! ! *** finish ! ! call Done( tmmd, status ) ! call Done( levi, status ) ! call Done( levi_ec, status ) ! call Done( lli ) ! ! ! -- end ! ! \ev ! ! ! !INTRODUCTION: Rcfile ! ! \bv ! ! ! ! ! Meteo files are linked to or unpacked in a buffer directory. ! ! o Set the clean flag (T|F) such that files that have not been accessed ! ! for a long time are removed if a maximum buffer usage is exceeded. ! ! o specify a maximum size in Mb ! ! ! tmm.dir : ${RUNDIR}/tmm-buf ! tmm.dir.clean : T ! tmm.dir.size : 500 ! ! ! ! ! TMM requires keys on how to form meteo for a certain region. ! ! A key should be defined for each region, names are in 'dims_grid.F90' ! ! For example: ! ! ! ! tmpp:od-fc-ml60-glb3x2 ! ! Read global 3x2, 60 level files produced by TMPP. ! ! Optionally, the meteo is combined over levels or grid cells. ! ! The files are expected to be present in the buffer directory ! ! specified below after 'tmm.buf.dir' . ! ! To have the appropriate files installed at the begin of a run, ! ! use the 'tmm.setup.*' stuff below. ! ! ! ! tmppS:od-fc-ml60-glb3x2 ! ! Idem, but also calls a script to search for an appropriate file ! ! from within the fortran program. ! ! The system call to this script turned out to be rather slow. ! ! This source type should be avoided therefore, but might be ! ! very usefull in case of limitted disk space. ! ! ! ! prism: ! ! Receive meteo from the prism coupler. ! ! ! tmm.sourcekey.glb6x4 : tmpp:od-fc-ml60-glb3x2 ! tmm.sourcekey.eur3x2 : tmpp:od-fc-ml60-glb3x2 ! tmm.sourcekey.eur1x1 : tmpp:od-fc-tropo25-eur1x1 ! ! ! ! ! Meteo files could be setup before the actual program is started. ! ! Fill the following settings: ! ! o set the apply flag apply this feature or not (T|F) ! ! o specify a list of meteo files to be installed (spm,uvsp, etc) ! ! o specify a list of meteo sources (od-fc-ml60-glb3x2 etc) ! ! o specify wether message are printend or not (T|F) ! ! ! tmm.setup.apply : T ! tmm.setup.files : spm uvsp t q cld sub surf ! tmm.setup.sources : od-fc-ml60-glb3x2 ! tmm.setup.verbose : T ! ! ! ! ! Archive(s) to be searched for monthly tar files. ! ! If more than one is specified (space seperated list), ! ! multiple directories are examined. ! ! o disk archives ! ! o tape archives ecfs/mos ('massive-storage-system') ! ! ! tmm.search.disk : ${DATADIR}/meteo ! tmm.search.mss : /nlh/TM/meteo ! ! \ev ! ! ! !INTRODUCTION: Source and scripts ! ! \bv ! ! tmm.f90 : Main routines and collecting data structure. ! Provides access to a list of open meteo files. ! ! tmm_mf.f90 : Search, open, close a meteo file. ! Calls shell script 'bin/tmm_getmeteo'. ! Calls specific routines for hdf/etc files. ! ! tmm_mf_hdf.f90 : Read fields from hdf files. ! ! tmm_param.f90 : Parameter specific stuff: ! what to do with temperature fields, ! what to do with mass fluxes etc. ! ! \ev ! !EOI ! ! ! The spm/spmid problem ... ! ! Glossary: ! spm : surface pressures for 00, 06, ...; ! read from 'spm_' files, computed from ecmwf lnsp ! spmid : average of surface pressures at begin/end interval ! ! To have the same algorithm as in TMPP, ! the following procedure should be used: ! ! 1) use ReadUVSP to get sp1, spm, sp2 ! spm is now in fact a spmid field (average of sp1 and sp2); ! in future it should be the real spm: ! ! call ReadUVSP( tmm, archivekey, tday, t1, t2, & ! lli, levi, sp1, spm, sp2, pu, pv, status ) ! ! 2) use this spm in calls to ReadField: ! ! call ReadField( tmm, archivekey, paramkey, tday, t1, t2, & ! lli, nuv, levi, spm, ll, status ) ! ! Examples of current implementation: ! ! 1) Temperature 6x4x25 from 3x2x60 ! read spm 3x2 from file ! horizontal mass average to 6x4x60 using spm 3x2 from ! vertical combination to 6x4x25 using provided sp 3x2 ! ! 2) Temperature 6x4x25 from N80x60 ! read spm N80 ! horizontal mass average to 6x4x60 using spm N80 and provided sp 3x2 ! --> should be spm 3x2 from spm N80 ! !############################################################################### ! #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 "tmm.inc" ! !############################################################################### module TMM use GO , only : gol, goPr, goErr, goBug, goLabel use GO , only : TDate, wrtgol use Grid , only : TllGridInfo, TggGridInfo, TshGridInfo, TLevelInfo use tmm_mf , only : TMeteoFile use tmm_info, only : TMeteoInfo, SetHistory, AddHistory implicit none ! --- in/out ------------------------------------- private public :: TTmMeteo public :: Init, Done public :: ReadField public :: Read_SP, Read_TQ public :: Read_Convec, Read_CloudCovers public :: Read_Diffus public :: Read_SR_OLS ! public :: ReadEqvLat #ifdef with_prism public :: TMM_READ_UVW #endif public :: WriteField public :: TMeteoInfo, SetHistory, AddHistory ! --- const --------------------------------------- character(len=*), parameter :: mname = 'tmm' ! maximum number of opened meteo files integer, parameter :: maxmf = 200 ! --- types --------------------------------------- type TTmMeteo ! rc file with paths etc character(len=256) :: rcfilename ! paths character(len=256) :: input_dir, output_dir ! list of meteo files integer :: nmf type(TMeteoFile) :: mf(maxmf) ! ! buffers for latest read raw field ! logical :: buf_filled ! character(len=10) :: buf_archivekey character(len=10) :: buf_paramkey type(TDate) :: buf_tday, buf_t1, buf_t2 character(len=1) :: buf_nuv character(len=1) :: buf_nw type(TMeteoInfo) :: buf_tmi ! character(len=2) :: buf_gridtype ! type(TllGridInfo) :: buf_lli real, pointer :: buf_ll(:,:,:) real, pointer :: buf_sp_ll(:,:) ! type(TggGridInfo) :: buf_ggi real, pointer :: buf_gg(:,:) real, pointer :: buf_sp_gg(:) ! type(TshGridInfo) :: buf_shi complex, pointer :: buf_sh(:,:) complex, pointer :: buf_lnsp_sh(:) ! type(TLevelInfo) :: buf_levi ! ! storage for horizontal interpol real, pointer :: llX(:,:,:) end type TTmMeteo ! --- interfaces ------------------------------------- interface Init module procedure tmm_Init end interface interface Done module procedure tmm_Done end interface interface SelectMF module procedure tmm_SelectMF end interface interface ReadField module procedure tmm_ReadField_2d module procedure tmm_ReadField_3d end interface ! interface ReadUVSP ! module procedure tmm_ReadUVSP ! end interface interface Read_SP module procedure tmm_Read_SP end interface interface Read_TQ module procedure tmm_Read_TQ end interface interface Read_Convec module procedure tmm_Read_Convec end interface interface Read_Diffus module procedure tmm_Read_Diffus end interface interface Read_CloudCovers module procedure tmm_Read_CloudCovers end interface interface Read_SR_OLS module procedure tmm_Read_SR_OLS end interface ! interface ReadEqvLat ! module procedure tmm_ReadEqvLat ! end interface interface WriteField module procedure tmm_WriteField_2d module procedure tmm_WriteField_3d end interface ! --- var -------------------------------------- ! timer id's: integer :: itim_fillbuffer integer :: itim_readfield_2d integer :: itim_readfield_3d integer :: itim_transform_2d integer :: itim_transform_3dh integer :: itim_transform_3dv contains ! =================================================================== subroutine tmm_Init( tmm, rcF, status ) use PArray , only : pa_Init use GO , only : NewDate use GO , only : TrcFile, ReadRc use GO , only : GO_Timer_Def #ifdef with_tmm_tm5 use tmm_mf_tm5_nc , only : TMM_MF_TM5_NC_Init #endif #ifdef with_prism use tmm_mf_prism, only : mfPrism_Init #endif ! --- in/out ---------------------------------- type(TTmMeteo), intent(inout) :: tmm type(TrcFile), intent(in) :: rcF integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Init' ! --- local ----------------------------------- ! --- begin ----------------------------------- ! store rc file name tmm%rcfilename = trim(rcf%fname) ! read paths call ReadRc( rcF, 'tmm.dir', tmm%input_dir, status ) IF_NOTOK_RETURN(status=1) ! read output path: set to a dummy value if key not found in rcfiles, ! since probably no meteo output requested anyway in this case ... call ReadRc( rcF, 'tmm.output.dir', tmm%output_dir, status, & default='/no/tmm/output/dir/specified/' ) ! no files open yet tmm%nmf = 0 ! buffer empty tmm%buf_filled = .false. tmm%buf_archivekey = 'none' tmm%buf_paramkey = 'none' tmm%buf_tday = NewDate(time5=(/0001,01,01,00,00/)) tmm%buf_t1 = NewDate(time5=(/0001,01,01,00,00/)) tmm%buf_t2 = NewDate(time5=(/9999,12,31,00,00/)) tmm%buf_nuv = 'n' tmm%buf_nw = 'n' tmm%buf_gridtype = 'xx' call pa_Init( tmm%buf_ll ) call pa_Init( tmm%buf_sp_ll ) call pa_Init( tmm%buf_gg ) call pa_Init( tmm%buf_sp_gg ) call pa_Init( tmm%buf_sh ) call pa_Init( tmm%buf_lnsp_sh ) ! init temp grid on destination grid but raw levels call pa_Init( tmm%llX ) #ifdef with_prism ! init prism stuff: call mfPrism_Init( status ) IF_NOTOK_RETURN(status=1) #endif #ifdef with_tmm_tm5 ! init input of TM5/NetCDF files: call TMM_MF_TM5_NC_Init( rcf, status ) IF_NOTOK_RETURN(status=1) #endif ! init timers: call GO_Timer_Def( itim_fillbuffer , 'tmm fill buffer', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_readfield_2d , 'tmm readfield 2D', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_readfield_3d , 'tmm readfield 3D', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_transform_2d , 'tmm transform 2D', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_transform_3dh, 'tmm transform 3D hor' , status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_transform_3dv, 'tmm transform 3D vert', status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine tmm_Init ! *** subroutine tmm_Done( tmm, status ) use PArray , only : pa_Done use tmm_mf , only : Opened, Done #ifdef with_tmm_tm5 use tmm_mf_tm5_nc , only : TMM_MF_TM5_NC_Done #endif #ifdef with_prism use tmm_mf_prism, only : mfPrism_Done #endif ! --- in/out ---------------------------------- type(TTmMeteo), intent(inout) :: tmm integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Done' ! --- local ---------------------------------------- integer :: imf ! --- begin ------------------------------------------- #ifdef with_tmm_tm5 ! init input of TM5/NetCDF files: call TMM_MF_TM5_NC_Done( status ) IF_NOTOK_RETURN(status=1) #endif #ifdef with_prism ! done with prism stuff: call mfPrism_Done( status ) IF_NOTOK_RETURN(status=1) #endif ! loop over all available meteo files do imf = 1, tmm%nmf ! in use ? if ( .not. Opened( tmm%mf(imf) ) ) cycle ! close ... call Done( tmm%mf(imf), status ) IF_NOTOK_RETURN(status=1) end do ! clear buffer call ClearBuffer( tmm, status ) IF_NOTOK_RETURN(status=1) call pa_Done( tmm%buf_ll ) call pa_Done( tmm%buf_sp_ll ) call pa_Done( tmm%buf_gg ) call pa_Done( tmm%buf_sp_gg ) call pa_Done( tmm%buf_sh ) call pa_Done( tmm%buf_lnsp_sh ) ! clear temp grid call pa_Done( tmm%llX ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ok status = 0 end subroutine tmm_Done ! *** ! Select meteo file for this param and time, ! search for a new file if necessary subroutine tmm_SelectMF( tmm, io, archivekey, paramkey, tday, t1, t2, imf, status ) use tmm_mf, only : Init, Done, Opened, CheckTime, CheckParam, SetupInput use tmm_mf, only : SetupInput, SetupOutput ! --- in/out ---------------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=1), intent(in) :: io character(len=*), intent(in) :: archivekey character(len=*), intent(in) :: paramkey type(TDate), intent(in) :: tday, t1, t2 integer, intent(out) :: imf integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_SelectMF' ! --- local ---------------------------------------- integer :: i ! --- begin ---------------------------------------- !write (gol,'(" DDD selectmf paramkey : ==",a,"==")')paramkey; call goBug !write (gol,*) 'DDD selectmf archivekey : ', trim(archivekey); call gobug !call wrtgol( ' DDD selectmf tday : ', tday ); call gobug !call wrtgol( ' DDD selectmf t1 - t2 : ', t1, ' - ', t2 ); call gobug ! not found yet ... imf = -1 ! loop over currently used meteo files do i = 1, tmm%nmf ! in use ? if ( .not. Opened( tmm%mf(i) ) ) cycle ! test on requested grid and param: call CheckParam( tmm%mf(i), io, archivekey, paramkey, status ) !write (gol,*) 'DDD selectmf chk param : ', i, trim(tmm%mf(i)%filename), status; call gobug if ( status == 0 ) then ! param included: leave imf = i exit else if ( status < 0 ) then ! param not included, try next cycle else ! error ... TRACEBACK; status=1; return end if end do ! time ok ? close if data in file is too old: if ( imf > 0 ) then ! time ok ? call CheckTime( tmm%mf(imf), t1, t2, status ) !write (gol,*) 'DDD selectmf chk time : ', status; call gobug if ( status == 0 ) then ! mf includes [t1,t2]; return with this imf status = 0; return else if ( status < 0 ) then ! mf does not include [t1,t2]; close file !write (gol,*) 'DDD selectmf close : ', imf, trim(tmm%mf(imf)%filename); call gobug call Done( tmm%mf(imf), status ) IF_NOTOK_RETURN(status=1) ! not in use anymore ... imf = -1 else TRACEBACK; status=1; return end if end if !write (gol,*) 'DDD selectmf imf : ', imf; call gobug ! open new meteo file ? if ( imf < 0 ) then ! search first available empty mf do i = 1, tmm%nmf if ( .not. Opened(tmm%mf(i)) ) then imf = i exit end if end do ! next number ? if ( imf < 0 ) then tmm%nmf = tmm%nmf + 1 if ( tmm%nmf > maxmf ) then write (gol,'("Tried to init meteo file beyond maximum number: ",i6)') maxmf; call goErr write (gol,'("Initialized files:")'); call goErr do i = 1, maxmf write (gol,'(" ",a)') trim(tmm%mf(i)%filename); call goErr end do write (gol,'("Please increase parameter `maxmf` in ",a)') mname; call goErr TRACEBACK; status=1; return end if imf = tmm%nmf end if ! start new mf ... call Init( tmm%mf(imf), io, status ) IF_NOTOK_RETURN(status=1) end if ! input or output ? select case ( io ) case ( 'i' ) ! input ! open file, store description, etc call SetupInput( tmm%mf(imf), archivekey, paramkey, tday, t1, t2, & tmm%rcfilename, tmm%input_dir, status ) IF_NOTOK_RETURN(status=1) case ( 'o' ) ! output ! open file, store description, etc call SetupOutput( tmm%mf(imf), archivekey, paramkey, tday, t1, t2, & tmm%rcfilename, tmm%output_dir, status ) IF_NOTOK_RETURN(status=1) case default write (gol,'("unsupported io `",a,"`")') io; call goErr TRACEBACK; status=1; return end select ! ok status = 0 end subroutine tmm_SelectMF ! ================================================================== ! === ! === buffer ! === ! ================================================================== subroutine FillBuffer( tmm, archivekey, paramkey, unit, tday, t1, t2, nuv, nw, status ) use GO , only : TDate, operator(==) use GO , only : GO_Timer_Start, GO_Timer_End use tmm_mf, only : ReadRecord ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey character(len=*), intent(in) :: paramkey character(len=*), intent(in) :: unit type(TDate), intent(in) :: tday, t1, t2 character(len=1), intent(in) :: nuv character(len=1), intent(in) :: nw integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/FillBuffer' ! --- local ------------------------------- integer :: imf ! --- begin ------------------------------- ! start timing call GO_Timer_Start( itim_fillbuffer, status ) IF_NOTOK_RETURN(status=1) ! is requested field already in buffer? if ( (archivekey == tmm%buf_archivekey) .and. (paramkey == tmm%buf_paramkey) .and. & (nuv == tmm%buf_nuv) .and. (nw == tmm%buf_nw) .and. & (tday == tmm%buf_tday) .and. (t1 == tmm%buf_t1) .and. (t2 == tmm%buf_t2) ) then call GO_Timer_End( itim_fillbuffer, status ) IF_NOTOK_RETURN(status=1) return end if ! select (eventually retrieve first) the meteo file with this param: call SelectMF( tmm, 'i', archivekey, paramkey, tday, t1, t2, imf, status ) IF_NOTOK_RETURN(status=1) ! clear buffer call ClearBuffer( tmm, status ) IF_NOTOK_RETURN(status=1) ! fill keys: tmm%buf_archivekey = archivekey tmm%buf_paramkey = paramkey tmm%buf_t1 = t1 tmm%buf_t2 = t2 tmm%buf_tday = tday tmm%buf_nuv = nuv tmm%buf_nw = nw ! read field: call ReadRecord( tmm%mf(imf), paramkey, unit, tday, t1, t2, nuv, nw, & tmm%buf_gridtype, tmm%buf_levi, & tmm%buf_lli, tmm%buf_ll, tmm%buf_sp_ll, & tmm%buf_ggi, tmm%buf_gg, tmm%buf_sp_gg, & tmm%buf_shi, tmm%buf_sh, tmm%buf_lnsp_sh, & tmm%buf_tmi, status ) IF_NOTOK_RETURN(status=1) ! some data present .. tmm%buf_filled = .true. ! end timing: call GO_Timer_End( itim_fillbuffer, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine FillBuffer ! *** subroutine ClearBuffer( tmm, status ) use PArray , only : pa_Done use GO , only : goErr use Grid , only : Done use tmm_info, only : Done ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/ClearBuffer' ! --- begin ------------------------------- if ( tmm%buf_filled ) then call Done( tmm%buf_tmi, status ) IF_NOTOK_RETURN(status=1) if ( tmm%buf_gridtype == 'll' ) then call Done( tmm%buf_lli, status ) IF_NOTOK_RETURN(status=1) call pa_Done( tmm%buf_ll ) call pa_Done( tmm%buf_sp_ll ) end if if ( tmm%buf_gridtype == 'gg' ) then call Done( tmm%buf_ggi, status ) IF_NOTOK_RETURN(status=1) call pa_Done( tmm%buf_gg ) call pa_Done( tmm%buf_sp_gg ) end if if ( tmm%buf_gridtype == 'sh' ) then call Done( tmm%buf_shi ) call pa_Done( tmm%buf_sh ) call pa_Done( tmm%buf_lnsp_sh ) end if call Done( tmm%buf_levi, status ) IF_NOTOK_RETURN(status=1) end if ! ok status = 0 end subroutine ClearBuffer ! *** subroutine CheckBuffer( tmm, status, gridtype, np, shT, nlev ) ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm integer, intent(out) :: status character(len=*), intent(in), optional :: gridtype integer, intent(in), optional :: np integer, intent(in), optional :: shT integer, intent(in), optional :: nlev ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/CheckBuffer' ! --- begin ------------------------------- if ( .not. tmm%buf_filled ) then write (gol,'("buffer not filled ...")'); call goErr TRACEBACK; status=1; return end if if ( present(gridtype) ) then if ( tmm%buf_gridtype /= gridtype ) then write (gol,'("unexpected gridtype in buffer:")'); call goErr write (gol,'(" buffer : ",a)') tmm%buf_gridtype; call goErr write (gol,'(" expected : ",a)') gridtype; call goErr TRACEBACK; status=1; return end if end if if ( present(nlev) ) then if ( tmm%buf_levi%nlev /= nlev ) then write (gol,'("unexpected number of levels in buffer:")'); call goErr write (gol,'(" buffer : ",i4)') tmm%buf_levi%nlev ; call goErr write (gol,'(" expected : ",i4)') nlev ; call goErr TRACEBACK; status=1; return end if end if if ( present(np) ) then if ( tmm%buf_ggi%np /= np ) then write (gol,'("unexpected grid size in buffer:")'); call goErr write (gol,'(" buffer ggi%np : ",i6)') tmm%buf_ggi%np; call goErr write (gol,'(" expected : ",i6)') np; call goErr TRACEBACK; status=1; return end if end if if ( present(shT) ) then if ( tmm%buf_shi%T /= shT ) then write (gol,'("unexpected grid size in buffer:")'); call goErr write (gol,'(" buffer shi%shT : ",i6)') tmm%buf_shi%T; call goErr write (gol,'(" expected : ",i6)') shT; call goErr TRACEBACK; status=1; return end if end if ! ok status = 0 end subroutine CheckBuffer ! *** subroutine Transform2D( tmm, lli, nuv, tmi, status ) ! use PArray , only : pa_SetShape use GO , only : GO_Timer_Start, GO_Timer_End use Grid , only : Tgg2llFracs, Init, Done use Grid , only : Interpol, FillGrid use Grid , only : FracSum use tmm_param , only : GetCombineKeys use tmm_info , only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm ! meteo buffer (incl. source data[in], target data[out], source grid) type(TllGridInfo), intent(in) :: lli ! grid of target data character(len=1), intent(in) :: nuv ! horizontal staggering of target data type(TMeteoInfo), intent(out) :: tmi ! integer, intent(out) :: status ! ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/Transform2D' ! --- local ------------------------------- character(len=10) :: hcomb, vcomb type(Tgg2llFracs) :: gg2ll ! --- begin ---------------------------------- ! start timing: call GO_Timer_Start( itim_transform_2d, status ) IF_NOTOK_RETURN(status=1) ! copy info: tmi = tmm%buf_tmi ! set shape of target grid: ! Note: switched from pa_SetShape to allocate when developping the parallel IO reading of nc met fields. if (associated (tmm%llX)) deallocate(tmm%llX) select case ( nuv ) case ( 'n' ) allocate( tmm%llX( lli%nlon , lli%nlat , 1 ) ) case ( 'u' ) allocate( tmm%llX( lli%nlon+1, lli%nlat , 1 ) ) case ( 'v' ) allocate( tmm%llX( lli%nlon , lli%nlat+1, 1 ) ) case default write (gol,'("unsupported nuv `",a,"`")') nuv TRACEBACK; status=1; return end select ! fill with zero's for safety: tmm%llX = 0.0 ! define how to combine horizontal cells and vertical levels: call GetCombineKeys( hcomb, vcomb, tmm%buf_paramkey, status ) IF_NOTOK_RETURN(status=1) ! transform raw field to ll : select case ( tmm%buf_gridtype ) ! ~~~ lat/lon ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'll' ) ! copy or combine horizontal cells from buffer into lli/llX : call FillGrid( lli, nuv, tmm%llX(:,:,1), & tmm%buf_lli, tmm%buf_nuv, tmm%buf_ll(:,:,1), & hcomb, status ) IF_ERROR_RETURN(status=1) ! Catch cases of incomplete mapping if (status==-1)then write(gol,*) " "; call goErr write(gol,*) "Only part of target array was filled!"; call goErr write(gol,*) "Make sure that the number of processors divides the number of grid boxes."; call goErr write(gol,*) " "; call goErr TRACEBACK; status=1; return endif ! ~~~ gg ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'gg' ) ! determine fraction call Init( gg2ll, tmm%buf_ggi, lli, status ) IF_NOTOK_RETURN(status=1) ! deceide how to interpolate given hcomb key: select case ( hcomb ) case ( 'area-aver' ) ! take fractions of overlapping cells call FracSum( gg2ll, tmm%buf_gg(:,1), tmm%llX(:,:,1), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) case default write (gol,'("unsupported horizonal combination key for gg :",a)') hcomb; call goErr write (gol,'("TIP: hcomb not set for this variable in module tmm_param ?")'); call goErr TRACEBACK; status=1; return end select ! done call Done( gg2ll ) ! ~~~ sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'sh' ) ! interpolate field in shi/sh into lli/ll : call Interpol( tmm%buf_shi, tmm%buf_sh(:,1), lli, tmm%llX(:,:,1), status ) IF_NOTOK_RETURN(status=1) case default write (gol,'("unsupported input grid type `",a,"`")') tmm%buf_gridtype; call goErr TRACEBACK; status=1; return end select ! fill history: call Init( tmi, tmm%buf_tmi, status ) call AddHistory( tmi, 'oper==hcomb,'//trim(hcomb), status ) ! end timing: call GO_Timer_End( itim_transform_2d, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine Transform2D ! *** subroutine Transform3Dh( tmm, lli, nuv, levi, nw, sp, tmi, status ) use PArray , only : pa_SetShape use GO , only : GO_Timer_Start, GO_Timer_End use Grid , only : FillLevels, FillGrid, AreaOper use Grid , only : Tgg2llFracs, Init, Done use Grid , only : IntArea use Grid , only : FracSum use tmm_param , only : GetCombineKeys use tmm_info , only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm type(TllGridInfo), intent(in) :: lli character(len=1), intent(in) :: nuv type(TLevelInfo), intent(in) :: levi character(len=1), intent(in) :: nw real, intent(out) :: sp(:,:) type(TMeteoInfo), intent(out) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/Transform3Dh' ! --- local ------------------------------- integer :: nlon, nlat, nlev character(len=10) :: hcomb, vcomb integer :: l type(Tgg2llFracs) :: gg2ll real, allocatable :: dp_ll(:,:) real, allocatable :: dp_llX(:,:) real, allocatable :: dp_gg(:) ! --- begin ---------------------------------- ! start timing: call GO_Timer_Start( itim_transform_3dh, status ) IF_NOTOK_RETURN(status=1) ! set shape of target grid: select case ( nuv ) case ( 'n' ) nlon = lli%nlon nlat = lli%nlat case ( 'u' ) nlon = lli%nlon+1 nlat = lli%nlat case ( 'v' ) nlon = lli%nlon nlat = lli%nlat+1 case default write (gol,'("unsupported nuv `",a,"`")') nuv; call goErr TRACEBACK; status=1; return end select select case ( nw ) case ( 'n' ) nlev = tmm%buf_levi%nlev case ( 'w' ) nlev = tmm%buf_levi%nlev+1 case default write (gol,'("unsupported nw `",a,"`")') nw; call goErr TRACEBACK; status=1; return end select call pa_SetShape( tmm%llX, nlon, nlat, nlev ) tmm%llX = 0.0 ! define how to combine horizontal cells and vertical levels: call GetCombineKeys( hcomb, vcomb, tmm%buf_paramkey, status ) IF_NOTOK_RETURN(status=1) ! transform raw field to ll : select case ( tmm%buf_gridtype ) ! ~~~ lat/lon ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'll' ) !! check size of buffer array !if ( size(tmm%buf_ll,3) /= nlev ) then ! write (gol,'("buffer array does not match with level definition:")') ! write (gol,'(" nw : ",a )') nw ! write (gol,'(" buf_levi nlev : ",i3)') tmm%buf_levi%nlev ! write (gol,'(" buf_ll levels : ",i3)') size(tmm%buf_ll,3) ! TRACEBACK; status=1; return !end if ! convective fluxes have only lowest layers ... if ( size(tmm%buf_ll,3) > nlev ) then write (gol,'("buffer array has more layers than implied by nw and level definition:")'); call goErr write (gol,'(" nw : ",a )') nw; call goErr write (gol,'(" buf_levi nlev : ",i3)') tmm%buf_levi%nlev; call goErr write (gol,'(" buf_ll levels : ",i3)') size(tmm%buf_ll,3); call goErr TRACEBACK; status=1; return end if ! surface pressure: aver horizontal cells from buffer into lli/llX : call FillGrid( lli, nuv, sp, & tmm%buf_lli, tmm%buf_nuv, tmm%buf_sp_ll, & 'area-aver', status ) IF_ERROR_RETURN(status=1) ! Catch cases of incomplete mapping if (status==-1)then write(gol,*) " "; call goErr write(gol,*) "Only part of target array was filled!"; call goErr write(gol,*) "Make sure that the number of processors divides the number of grid boxes."; call goErr write(gol,*) " "; call goErr TRACEBACK; status=1; return endif ! deceide how to interpolate given hcomb key; ! for most keys, let 'FillGrid' determine what to do ... select case ( hcomb ) case ( 'mass-aver' ) ! check ... if ( nw /= 'n' ) then write (gol,'("nw should be `n` for mass aver ...")'); call goErr TRACEBACK; status=1; return end if ! ! mass weighted fractions: ! ! sum_i f_i dp_i/g dA_i ! ---------------------- ! sum_i dp_i/g dA_i ! ! temporary pressure field: allocate( dp_ll(tmm%buf_lli%nlon,tmm%buf_lli%nlat) ) ! loop over layers !do l = 1, nlev do l = 1, size(tmm%buf_ll,3) ! pressure gradients in this layer: dp_ll = tmm%buf_levi%da(l) + tmm%buf_levi%db(l) * tmm%buf_sp_ll call AreaOper( tmm%buf_lli, dp_ll, '*', 'm2', status ) IF_NOTOK_RETURN(status=1) ! copy or combine horizontal cells from buffer into lli/llX; ! combining cells is weighted by 'mass' dp_ll : call FillGrid( lli, nuv, tmm%llX(:,:,l), & tmm%buf_lli, tmm%buf_nuv, tmm%buf_ll(:,:,l), & 'weight', status, dp_ll ) IF_ERROR_RETURN(status=1) ! Catch cases of incomplete mapping if (status==-1)then write(gol,*) " "; call goErr write(gol,*) "Only part of target array was filled!"; call goErr write(gol,*) "Make sure that the number of processors divides the number of grid boxes."; call goErr write(gol,*) " "; call goErr TRACEBACK; status=1; return endif end do ! clear: deallocate( dp_ll ) case default ! loop over layers in ll array: !do l = 1, nlev do l = 1, size(tmm%buf_ll,3) ! copy or combine horizontal cells from buffer into lli/llX; ! pass hcomb to FillGrid to define how cells are combined: call FillGrid( lli, nuv, tmm%llX(:,:,l), & tmm%buf_lli, tmm%buf_nuv, tmm%buf_ll(:,:,l), & hcomb, status ) IF_ERROR_RETURN(status=1) ! Catch cases of incomplete mapping if (status==-1)then write(gol,*) " "; call goErr write(gol,*) "Only part of target array was filled!"; call goErr write(gol,*) "Make sure that the number of processors divides the number of grid boxes."; call goErr write(gol,*) " "; call goErr TRACEBACK; status=1; return endif end do end select ! ~~~ gg ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'gg' ) ! check size of buffer array if ( size(tmm%buf_gg,2) /= nlev ) then write (gol,'("buffer array does not match with level definition:")'); call goErr write (gol,'(" nw : ",a )') nw; call goErr write (gol,'(" buf_levi nlev : ",i3)') tmm%buf_levi%nlev; call goErr write (gol,'(" buf_gg levels : ",i3)') size(tmm%buf_gg,2); call goErr TRACEBACK; status=1; return end if ! determine fraction call Init( gg2ll, tmm%buf_ggi, lli, status ) IF_NOTOK_RETURN(status=1) ! surface pressure: take fractions of overlapping cells call FracSum( gg2ll, tmm%buf_sp_gg, sp, status, 'area-aver' ) IF_NOTOK_RETURN(status=1) ! deceide how to interpolate given hcomb key: select case ( hcomb ) case ( 'area-aver' ) ! ! area weighted fractions: ! ! sum_i f_i dA_i ! --------------- ! sum_i dA_i ! ! loop over layers do l = 1, nlev ! take fractions of overlapping cells call FracSum( gg2ll, tmm%buf_gg(:,l), tmm%llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do case ( 'mass-aver' ) ! check ... if ( nw /= 'n' ) then write (gol,'("nw should be `n` for mass aver ...")'); call goErr TRACEBACK; status=1; return end if ! ! mass weighted fractions: ! ! sum_i f_i dp_i/g dA_i ! ---------------------- ! sum_i dp_i/g dA_i ! ! temporary pressure field: allocate( dp_gg(size(tmm%buf_sp_gg)) ) allocate( dp_llX(size(sp,1),size(sp,2)) ) ! loop over layers do l = 1, nlev ! pressure gradients in this layer: dp_gg = tmm%buf_levi%da(l) + tmm%buf_levi%db(l) * tmm%buf_sp_gg dp_llX = tmm%buf_levi%da(l) + tmm%buf_levi%db(l) * sp ! take fractions of overlapping cells: ! sum_i f_i dp_i dA_i call FracSum( gg2ll, tmm%buf_gg(:,l)*dp_gg, tmm%llX(:,:,l), status, 'area-sum' ) IF_NOTOK_RETURN(status=1) ! weight by total dp dA tmm%llX(:,:,l) = tmm%llX(:,:,l) / dp_llX call AreaOper( lli, tmm%llX(:,:,l), '/', 'm2', status ) IF_NOTOK_RETURN(status=1) end do ! clear: deallocate( dp_gg ) deallocate( dp_llX ) case default write (gol,'("unsupported horizonal combination key for gg :",a)') hcomb; call goErr TRACEBACK; status=1; return end select ! done call Done( gg2ll ) ! ~~~ sh ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'sh' ) ! check size of buffer array if ( size(tmm%buf_sh,2) /= nlev ) then write (gol,'("buffer array does not match with level definition:")'); call goErr write (gol,'(" nw : ",a )') nw; call goErr write (gol,'(" buf_levi nlev : ",i3)') tmm%buf_levi%nlev; call goErr write (gol,'(" buf_sh levels : ",i3)') size(tmm%buf_sh,2); call goErr TRACEBACK; status=1; return end if ! surface pressure: interpolate lnsp to target horizontal resolution: call IntArea( 'exp,aver', tmm%buf_shi, tmm%buf_lnsp_sh, lli, sp, status ) IF_NOTOK_RETURN(status=1) ! deceide how to interpolate given hcomb key: select case ( hcomb ) !case ( 'area-aver' ) ! ! ! area average over sepectral fields: ! call IntArea( 'aver', tmm%buf_shi, tmm%buf_sh, lli, tmm%llX, status ) ! IF_NOTOK_RETURN(status=1) case ( 'mass-aver' ) ! check ... if ( nw /= 'n' ) then write (gol,'("nw should be `n` for mass aver ...")'); call goErr TRACEBACK; status=1; return end if ! mass average over sepectral fields: call IntArea( '[F*(da+db*exp(H))*cos]/[()*cos]', & tmm%buf_shi, nlev, tmm%buf_sh, & tmm%buf_lnsp_sh, tmm%buf_levi%da, tmm%buf_levi%db, & lli, tmm%llX, status ) IF_NOTOK_RETURN(status=1) case default write (gol,'("unsupported horizonal combination key for sh :",a)') hcomb; call goErr TRACEBACK; status=1; return end select ! ~~~ ?? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported input grid type `",a,"`")') tmm%buf_gridtype; call goErr TRACEBACK; status=1; return end select ! fill history: call Init( tmi, tmm%buf_tmi, status ) call AddHistory( tmi, 'oper==hcomb,'//trim(hcomb), status ) ! end timing: call GO_Timer_End( itim_transform_3dh, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine Transform3Dh ! *** subroutine Transform3Dv( tmm, levi, nw, sp, ll, tmi, status ) use PArray , only : pa_SetShape use GO , only : GO_Timer_Start, GO_Timer_End use Grid , only : FillLevels, FillGrid, AreaOper use Grid , only : Tgg2llFracs, Init, Done use Grid , only : IntArea use Grid , only : FracSum use tmm_param , only : GetCombineKeys use tmm_info , only : TMeteoInfo, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm type(TLevelInfo), intent(in) :: levi character(len=1), intent(in) :: nw real, intent(in) :: sp(:,:) real, intent(out) :: ll(:,:,:) type(TMeteoInfo), intent(inout) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/Transform3Dv' ! --- local ------------------------------- character(len=10) :: hcomb, vcomb ! --- begin ---------------------------------- ! start timing: call GO_Timer_Start( itim_transform_3dv, status ) IF_NOTOK_RETURN(status=1) ! define how to combine horizontal cells and vertical levels: call GetCombineKeys( hcomb, vcomb, tmm%buf_paramkey, status ) IF_NOTOK_RETURN(status=1) ! collect or copy levels, eventually reversed, from leviX/llX into levi/ll : !write (gol,'(a,": vertical ...")') rname call FillLevels( levi, nw, sp, ll, tmm%buf_levi, tmm%llX, vcomb, status ) IF_NOTOK_RETURN(status=1) ! add history: call AddHistory( tmi, 'oper==vcomb,'//trim(vcomb), status ) ! end timing: call GO_Timer_End( itim_transform_3dv, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine Transform3Dv ! ================================================================== ! === ! === read ll field ! === ! ================================================================== subroutine tmm_ReadField_2d( tmm, archivekey, paramkey, unit, tday, t1, t2, & lli, nuv, ll, tmi, status ) use PArray , only : pa_Init, pa_Done use GO , only : GO_Timer_Start, GO_Timer_End use tmm_info, only : TMeteoInfo ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm ! meteo buffer (all about the SRC, TGT data) character(len=*), intent(in) :: archivekey character(len=*), intent(in) :: paramkey character(len=*), intent(in) :: unit type(TDate), intent(in) :: tday, t1, t2 ! date, time interval of request type(TllGridInfo), intent(in) :: lli ! grid info of requested data. Tile only if //-IO. character(len=1), intent(in) :: nuv ! horizontal staggering of requested data real, intent(out) :: ll(:,:) ! read (and regridded) data. Bounds are lost. type(TMeteoInfo), intent(out) :: tmi ! info integer, intent(out) :: status ! return code ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_ReadField_2d' ! --- local ------------------------------- character(len=10) :: hcomb, vcomb character(len=32) :: paramkey_in ! --- begin ---------------------------------- call goLabel(rname) ! start timing: call GO_Timer_Start( itim_readfield_2d, status ) IF_NOTOK_RETURN(status=1) !DBG call wrtgol( ' tmm read 2D : '//trim(paramkey)//' (', tday, ') ', t1, ' - ', t2 ); call goPr ! check shape of target grid: if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) ) then write (gol,'("target array does not match with grid definition:")'); call goErr write (gol,'(" param : ",a )') paramkey; call goErr write (gol,'(" lli : ",i3," x ",i3)') lli%nlon, lli%nlat; call goErr write (gol,'(" nuv : ",a )') nuv; call goErr write (gol,'(" ll : ",i3," x ",i3)') shape(ll); call goErr TRACEBACK; status=1; return end if ! ! convert grid ! ! standard values? if ( trim(archivekey) == 'standard' ) then ! dummy value ll = 0.0 else ! translate? #ifdef with_ec_aver select case ( trim(paramkey) ) case ( 'lsp' ) ; paramkey_in = 'mlspr' case ( 'cp' ) ; paramkey_in = 'mcpr' case ( 'sf' ) ; paramkey_in = 'msr' case ( 'sshf' ) ; paramkey_in = 'msshf' case ( 'slhf' ) ; paramkey_in = 'mslhf' case ( 'ssrd' ) ; paramkey_in = 'msdwswrf' case ( 'strd' ) ; paramkey_in = 'msdwlwrf' case ( 'ssr' ) ; paramkey_in = 'msnswrf' case ( 'str' ) ; paramkey_in = 'msnlwrf' case ( 'ewss' ) ; paramkey_in = 'metss' case ( 'nsss' ) ; paramkey_in = 'mntss' case default paramkey_in = trim(paramkey) end select #else paramkey_in = trim(paramkey) #endif ! read raw field in buffer call FillBuffer( tmm, archivekey, paramkey_in, unit, tday, t1, t2, nuv, 'n', status ) IF_NOTOK_RETURN(status=1) ! horizontal interpolation: call Transform2D( tmm, lli, nuv, tmi, status ) IF_NOTOK_RETURN(status=1) ! expecting single level ... if ( size(tmm%llX,3) /= 1 ) then write (gol,'("expecting single level:")'); call goErr if ( trim(paramkey_in) /= trim(paramkey) ) then write (gol,'(" paramkey : ",a," (",a,")")') trim(paramkey_in), trim(paramkey); call goErr else write (gol,'(" paramkey : ",a)') trim(paramkey_in); call goErr end if write (gol,'(" levels : ",a)') size(tmm%llX,3); call goErr TRACEBACK; status=1; return end if ! extract first level ll = tmm%llX(:,:,1) end if ! ! unit conversions, truncations, etc ! select case ( paramkey ) case ( 'oro', 'cp', 'lsp' ) ! set minium; some negative values if made from spectral field: ll = max( 0.0, ll ) end select ! ! done ! ! end timing: call GO_Timer_End( itim_readfield_2d, status ) IF_NOTOK_RETURN(status=1) ! ok call goLabel() status = 0 end subroutine tmm_ReadField_2d ! *** subroutine tmm_ReadField_3d( tmm, archivekey, paramkey, unit, tday, t1, t2, & lli, nuv, levi, nw, sp, ll, tmi, status ) use GO , only : GO_Timer_Start, GO_Timer_End use Binas , only : p_global use tmm_info, only : TMeteoInfo ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey character(len=*), intent(in) :: paramkey, unit type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli character(len=1), intent(in) :: nuv type(TLevelInfo), intent(in) :: levi character(len=1), intent(in) :: nw real, intent(out) :: sp(:,:) real, intent(out) :: ll(:,:,:) type(TMeteoInfo), intent(out) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_ReadField_3d' ! --- begin ---------------------------------- call goLabel(rname) ! start timing: call GO_Timer_Start( itim_readfield_3d, status ) IF_NOTOK_RETURN(status=1) !DBG call wrtgol( ' tmm read3D : '//trim(paramkey)//' ', t1, ' - ', t2 ); call goPr ! check shape of target grid: if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) .or. & ((nuv == 'n') .and. ((size(sp,1) /= lli%nlon ) .or. (size(sp,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(sp,1) /= lli%nlon+1) .or. (size(sp,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(sp,1) /= lli%nlon ) .or. (size(sp,2) /= lli%nlat+1))) .or. & ((nw == 'n') .and. (size(ll,3) /= levi%nlev )) .or. & ((nw == 'w') .and. (size(ll,3) /= levi%nlev+1)) ) then write (gol,'("target arrays do not match with grid definition:")'); call goErr write (gol,'(" param : ",a )') paramkey ; call goErr write (gol,'(" lli : ",i3," x ",i3 )') lli%nlon, lli%nlat; call goErr write (gol,'(" nuv : ",a )') nuv; call goErr write (gol,'(" levi : ",i3 )') levi%nlev; call goErr write (gol,'(" nw : ",a )') nw; call goErr write (gol,'(" sp : ",i3," x ",i3 )') shape(sp); call goErr write (gol,'(" ll : ",i3," x ",i3," x ",i3)') shape(ll); call goErr TRACEBACK; status=1; return end if ! ! read, convert ! ! standard values? if ( trim(archivekey) == 'standard' ) then ! dummy values sp = p_global ll = 0.0 else ! read raw field in buffer call FillBuffer( tmm, archivekey, paramkey, unit, tday, t1, t2, nuv, nw, status ) IF_NOTOK_RETURN(status=1) ! 3d horizontal conversion: call Transform3Dh( tmm, lli, nuv, levi, nw, sp, tmi, status ) IF_NOTOK_RETURN(status=1) ! 3d vertical conversion: call Transform3Dv( tmm, levi, nw, sp, ll, tmi, status ) IF_NOTOK_RETURN(status=1) end if ! ! unit conversion, extreme values ! select case ( paramkey ) case ( 'Q', 'CLWC', 'CIWC' ) ! set minimum, stored values could be slightly negative ll = max( 0.0, ll ) ll = max( 0.0, ll ) end select ! ! done ! ! start timing: call GO_Timer_End( itim_readfield_3d, status ) IF_NOTOK_RETURN(status=1) ! ok call goLabel() status = 0 end subroutine tmm_ReadField_3d ! ================================================================== ! === ! === VO/D/LNSP -> pu/pv/sp ! === ! ================================================================== subroutine tmm_Read_SP( tmm, archivekey, paramname, paramunit, tday, t1, t2, & lli, sp, tmi, status ) use binas , only : p_global use GO , only : goSplitLine, Get use Grid , only : TshGrid, Init, Done, Set use Grid , only : IntArea use Grid , only : Tgg2llFracs, FracSum use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey, paramname, paramunit type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli real, intent(out) :: sp(:,:) ! Pa type(TMeteoInfo), intent(out) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_SP' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: hour type(Tgg2llFracs) :: gg2ll ! --- begin ------------------------------- call goLabel(rname) ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! standard ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'standard' ) !write (gol,'(" tmm read : sp standard global mean pressure")'); call goPr ! fill field with global mean pressure: sp = p_global ! set history info call Init( tmi, 'sp', 'Pa', status ) call AddHistory( tmi, 'standard', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from tmpp hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tmpp' ) ! select parameter given hour: call Get( t1, hour=hour ) if ( modulo(hour,6) == 3 ) then ! hour = 21/03/.. : uvsp files call ReadField( tmm, archivekey, 'sp', 'Pa', tday, t1, t1, & lli, 'n', sp, tmi, status ) IF_NOTOK_RETURN(status=1) else ! hour = 00/06/.. : spm files call ReadField( tmm, archivekey, 'spm', 'Pa', tday, t1, t1, & lli, 'n', sp, tmi, status ) IF_NOTOK_RETURN(status=1) end if ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from tm5 hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tm5-hdf', 'tm5-nc' ) ! sp always storred in 'sp' files ... call ReadField( tmm, archivekey, paramname, paramunit, tday, t1, t1, & lli, 'n', sp, tmi, status ) IF_NOTOK_RETURN(status=1) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! convert spectral lnsp ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tmpp', 'ncep-cdc', 'ncep-gfs', 'msc-data', 'prism' ) !DBG call wrtgol( ' tmm read : sp ', t1, ' - ', t2 ); call goPr ! read LNSP call FillBuffer( tmm, archivekey, 'LNSP', '1', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) !write (gol,*) " read LNSP and convert to SP!"; call goPr ! should be spectral ... if ( tmm%buf_gridtype /= 'sh' ) then write (gol,'("expecting sh field, not ",a)') tmm%buf_gridtype TRACEBACK; status=1; return end if ! aera average exp(lnsp) call IntArea( 'exp,aver', tmm%buf_shi, tmm%buf_sh(:,1), lli, sp, status ) IF_NOTOK_RETURN(status=1) ! set history info call Init( tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) call AddHistory( tmi, 'oper==exp,aver', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! convert gaussian sp or spectral lnsp ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tm5' ) select case ( paramname ) case ( 'sp', 'SP' ) !DBG call wrtgol( ' tmm read : sp ', t1, ' - ', t2 ); call goPr ! read LNSP call FillBuffer( tmm, archivekey, 'LNSP', '1', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! should be spectral ... if ( tmm%buf_gridtype /= 'sh' ) then write (gol,'("expecting sh field, not ",a)') tmm%buf_gridtype TRACEBACK; status=1; return end if ! aera average exp(lnsp) call IntArea( 'exp,aver', tmm%buf_shi, tmm%buf_sh(:,1), lli, sp, status ) IF_NOTOK_RETURN(status=1) ! set history info call Init( tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) call AddHistory( tmi, 'oper==exp,aver', status ) case ( 'sps', 'SPS' ) !DBG call wrtgol( ' tmm read : sps ', t1, ' - ', t2 ); call goPr ! read gg field call FillBuffer( tmm, archivekey, 'sp', 'Pa', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! should be gg ... if ( tmm%buf_gridtype /= 'gg' ) then write (gol,'("expecting gg field, not ",a)') tmm%buf_gridtype TRACEBACK; status=1; return end if ! determine fraction call Init( gg2ll, tmm%buf_ggi, lli, status ) IF_NOTOK_RETURN(status=1) ! take fractions of overlapping cells call FracSum( gg2ll, tmm%buf_gg(:,1), sp, status, 'area-aver' ) IF_NOTOK_RETURN(status=1) ! done call Done( gg2ll ) ! set history info call Init( tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) call AddHistory( tmi, 'oper==area-aver', status ) case default write (gol,'("unsupported param `",a,"` for source type `",a,"`")') & trim(paramname), trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! error ... ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ok call goLabel() status = 0 end subroutine tmm_Read_SP ! ***************************************************************** #ifdef with_prism SUBROUTINE TMM_READ_UVW( tmm, archivekey, tday, t1, t2, & lli, levi, & spu, mfu, mfu_tmi, & spv, mfv, mfv_tmi, & sp, mfw, tsp, tmi, & status ) use Binas , only : R => ae use Binas , only : g => grav use Binas , only : p_global use GO , only : goSplitLine use Grid , only : TshGrid, Init, Done, Set use Grid , only : IntLat, IntLon, vod2uv use Grid , only : FillLevels, Nabla, IntArea, AreaOper use tmm_info, only : TMeteoInfo, Init, Done, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(out) :: spu(:,:) , spv(:,:) ! Pa real, intent(out) :: mfu(:,:,:), mfv(:,:,:) ! kg/s ==> lower bounds become 1 !!!! type(TMeteoInfo), intent(out) :: mfu_tmi, mfv_tmi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: mfw(:,:,:) ! kg/s real, intent(out) :: tsp(:,:) ! tendency of surface pressure Pa/s type(TMeteoInfo), intent(out) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_UVW' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename type(TshGridInfo) :: shi integer :: shT integer :: nlev ! spectral fields type(TshGrid) :: LNSP_sh complex , pointer :: D_sh(:,:), VO_sh(:,:) complex , pointer :: U_sh(:,:), V_sh(:,:) complex , pointer :: Help_LNSP_sh(:,:) ! nabla.lnps type(TshGrid) :: NabLNSP_sh(2) ! integrated Omega arrays: real, pointer :: IIOmega (:,:,:) real, pointer :: IIOmega2(:,:,:) ! mfu/mfv/mfw on source levels real, allocatable :: mfuX(:,:,:) real, allocatable :: mfvX(:,:,:) real, pointer :: mfwX(:,:,:) ! loops etc integer :: l ! extra info type(TMeteoInfo) :: LNSP_tmi type(TMeteoInfo) :: vo_tmi, div_tmi, u_tmi, v_tmi ! temporary arrays real, allocatable :: tmp_sp(:,:,:) ! --- begin ------------------------------- call goLabel(rname) ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) case ( 'prism' ) !DBG call wrtgol( ' tmm read : mfuvw ', t1, ' - ', t2 ); call goPr ! read LNSP call FillBuffer( tmm, archivekey, 'LNSP', 'Pa', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh' ) IF_NOTOK_RETURN(status=1) ! extract grid size call Init( shi, tmm%buf_shi, status ) IF_NOTOK_RETURN(status=1) shT = tmm%buf_shi%T ! copy 1d spectral lnsp from buffer: call Init( LNSP_sh ) call Set( LNSP_sh, tmm%buf_shi%T, tmm%buf_sh(:,1) ) LNSP_tmi = tmm%buf_tmi ! read VO call FillBuffer( tmm, archivekey, 'VO', '1/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T ) IF_NOTOK_RETURN(status=1) ! extract grid size nlev = tmm%buf_levi%nlev ! copy 3d spectral field from buffer: allocate( VO_sh(shi%np,nlev) ) VO_sh = tmm%buf_sh vo_tmi = tmm%buf_tmi ! read D call FillBuffer( tmm, archivekey, 'D', '1/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d spectral field from buffer: allocate( D_sh(shi%np,nlev) ) D_sh = tmm%buf_sh div_tmi = tmm%buf_tmi ! ! compute U/V from VO/D ! allocate( U_sh(shi%np,nlev) ) allocate( V_sh(shi%np,nlev) ) allocate( Help_LNSP_sh(shi%np,1) ) !$OMP PARALLEL & !$OMP default ( none ) & !$OMP shared ( nlev, shi, VO_sh, D_sh, U_sh, V_sh ) & !$OMP private ( l ) do l = 1, nlev call vod2uv( shi, VO_sh(:,l), D_sh(:,l), shi, U_sh(:,l), V_sh(:,l) ) enddo !$OMP END PARALLEL ! history ... call Init( u_tmi, 'U', 'm/s', status, (/vo_tmi,div_tmi/) ) call Init( v_tmi, 'V', 'm/s', status, (/vo_tmi,div_tmi/) ) ! ! MFU = R/g int U (da+db*exp(LNSP)) / cos(lat) dlat ! allocate( mfuX(0:lli%nlon,lli%nlat,nlev) ) allocate( tmp_sp(0:lli%nlon,lli%nlat,1) ) ! integral over boundary: call IntLat( '(da+exp*db)/cos', shi, nlev, U_sh, LNSP_sh%c, & tmm%buf_levi%da, tmm%buf_levi%db, lli, mfuX, status ) IF_NOTOK_RETURN(status=1) mfuX = mfuX * R/g ! average surface pressure on boundary: Help_LNSP_sh(:,1)=LNSP_sh%c(:) call IntLat( 'exp(H),aver',shi,1, Help_LNSP_sh, LNSP_sh%c, (/0.0/), (/0.0/), lli, tmp_sp, status ) IF_NOTOK_RETURN(status=1) spu = tmp_sp(:,:,1) ! collect levels: call FillLevels( levi, 'n', spu, mfu, tmm%buf_levi, mfuX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! clear: deallocate( mfuX ) deallocate( tmp_sp ) ! info ... call Init( mfu_tmi, 'mfu', 'kg/s', status, (/u_tmi,LNSP_tmi/) ) call AddHistory( mfu_tmi, 'oper==intlat', status ) call AddHistory( mfu_tmi, 'oper==collectlevels', status ) ! ! MFV = R/g int V (da+db*exp(LNSP)) dlon ! allocate( mfvX(lli%nlon,0:lli%nlat,nlev) ) allocate( tmp_sp(lli%nlon,0:lli%nlat,1) ) ! integral over boundary: call IntLon( '(da+exp*db)',shi,nlev, V_sh, LNSP_sh%c, & tmm%buf_levi%da, tmm%buf_levi%db, lli, mfvX, status ) IF_NOTOK_RETURN(status=1) mfvX = mfvX * R/g ! average surface pressure on boundary: Help_LNSP_sh(:,1)=LNSP_sh%c(:) call IntLon( 'exp(H),aver', shi,1,Help_LNSP_sh, LNSP_sh%c, (/0.0/), (/0.0/), lli, tmp_sp, status ) IF_NOTOK_RETURN(status=1) spv = tmp_sp(:,:,1) ! collect levels: call FillLevels( levi, 'n', spv, mfv, tmm%buf_levi, mfvX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! clear: deallocate( Help_LNSP_sh) deallocate( mfvX ) deallocate( tmp_sp ) ! info ... call Init( mfv_tmi, 'mfv', 'kg/s', status, (/v_tmi,LNSP_tmi/) ) call AddHistory( mfv_tmi, 'oper==intlon', status ) call AddHistory( mfv_tmi, 'oper==collectlevels', status ) ! ! MFW ! allocate( IIOmega (lli%nlon,lli%nlat,nlev) ) allocate( IIOmega2(lli%nlon,lli%nlat,nlev) ) allocate( mfwX(lli%nlon,lli%nlat,nlev+1) ) ! ! int int D (da+db*exp(LNSP)) cos(lat) dA ! IIOmega = 0.0 call IntArea( 'F*(da+db*exp(H))*cos', shi, nlev, D_sh, LNSP_sh%c, LNSP_sh%c, & tmm%buf_levi%da, tmm%buf_levi%db, lli, IIOmega, status ) IF_NOTOK_RETURN(status=1) ! int int U dLNSP1 exp(LNSP) db / cos(lat) dA ! int int V dLNSP2 exp(LNSP) db / cos(lat) dA ! call Init( NabLNSP_sh(1) ) call Init( NabLNSP_sh(2) ) call Nabla( LNSP_sh, NabLNSP_sh ) ! compute nabla.lnsp ! integral over cell area; add contributions IIOmega2 = 0.0 call IntArea( 'F*G*(db*exp(H))/cos', shi, nlev, U_sh, NabLNSP_sh(1)%c, LNSP_sh%c, & tmm%buf_levi%da, tmm%buf_levi%db, lli, IIOmega2, status ) IF_NOTOK_RETURN(status=1) call IntArea( 'F*G*(db*exp(H))/cos', shi, nlev, V_sh, NabLNSP_sh(2)%c, LNSP_sh%c, & tmm%buf_levi%da, tmm%buf_levi%db, lli, IIOmega2, status ) IF_NOTOK_RETURN(status=1) ! deallocate call Done( NabLNSP_sh(1) ) call Done( NabLNSP_sh(2) ) ! ! column integrated Omega ! ! parent levels downwards or upwards ? if ( tmm%buf_levi%updo == 'd' ) then ! loop from top to bottom: do l = 1, nlev ! replace with contribution of current level: IIOmega(:,:,l) = (R**2)*IIOmega(:,:,l)/g + R*IIOmega2(:,:,l)/g ! add contribution of previous levels: if ( l > 1 ) then IIOmega(:,:,l) = IIOmega(:,:,l) + IIOmega(:,:,l-1) end if end do else ! loop from top to bottom: do l = nlev, 1, -1 ! replace with contribution of current level: IIOmega(:,:,l) = (R**2)*IIOmega(:,:,l)/g + R*IIOmega2(:,:,l)/g ! add contribution of levels above: if ( l < nlev ) then IIOmega(:,:,l) = IIOmega(:,:,l) + IIOmega(:,:,l+1) end if end do end if ! ! tendency of surface pressure ! ! dps 1 dp ! --- = - int nabla ( v ---- ) deta = - IIOmega(:,:,bot) ! dt eta=0 deta ! parent levels downwards or upwards ? if ( tmm%buf_levi%updo == 'd' ) then tsp = -1.0 * IIOmega(:,:,nlev) ! kg/s else tsp = -1.0 * IIOmega(:,:,1) ! kg/s end if ! convert to Pa/s : call AreaOper( lli, tsp, '/', 'm2', status ) ! kg/m2/s IF_NOTOK_RETURN(status=1) tsp = tsp * g ! Pa/s ! ! compute vertical flux: ! ! parent levels downwards or upwards ? if ( tmm%buf_levi%updo == 'd' ) then ! top hlev: mfwX(:,:,1) = 0.0 ! kg/s ! loop from top to bottom layer: do l = 1, nlev ! lay l bot hlev surflay bot hlev lay l bot hlev lay l bot hlev ! 2 .. nlev+1 nlev 0 .. nlev-1 1 .. nlev mfwX(:,:,l+1) = IIOmega(:,:,nlev) * tmm%buf_levi%b(l) - IIOmega(:,:,l) end do else ! top hlev: mfwX(:,:,nlev+1) = 0.0 ! kg/s ! loop from top to bottom layer: do l = nlev, 1, -1 ! lay l bot hlev surflay bot hlev lay l bot hlev lay l bot hlev ! 1 .. nlev nlev 0 .. nlev-1 1 .. nlev mfwX(:,:,l) = IIOmega(:,:,1) * tmm%buf_levi%b(l-1) - IIOmega(:,:,l) end do end if ! check: fluxh through bottom should be zero ... if ( maxval(abs(mfwX(:,:,1))) > 1.0 ) then write (gol,'("vert.flux through bottom half level should be zero ...")'); call goErr write (gol,'(" max value : ",es12.4)') maxval(abs(mfwX(:,:,1))); call goErr TRACEBACK; status=1; return end if ! create history call Init( tmi, 'mfw', 'kg/s', status, (/lnsp_tmi,div_tmi,u_tmi,v_tmi/) ) ! ! collect levels ! ! average surface pressure call IntArea( 'exp,aver', tmm%buf_shi, LNSP_sh%c, lli, sp, status ) IF_NOTOK_RETURN(status=1) ! combine levels etc call FillLevels( levi, 'w', sp, mfw, tmm%buf_levi, mfwX, 'bottom', status ) IF_NOTOK_RETURN(status=1) call AddHistory( tmi, 'oper==collectlevels', status ) ! ! upward flux ! ! flux should be upwards (in direction of increasing level): mfw = - mfw call AddHistory( tmi, 'oper==upwards', status ) ! ! done ! call Done( LNSP_sh ) deallocate( VO_sh ) deallocate( D_sh ) deallocate( U_sh ) deallocate( V_sh ) deallocate( IIOmega ) deallocate( IIOmega2 ) deallocate( mfwX ) call Done( lnsp_tmi, status ) call Done( vo_tmi, status ) call Done( div_tmi, status ) call Done( u_tmi, status ) call Done( v_tmi, status ) case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select call goLabel() status = 0 END SUBROUTINE TMM_READ_UVW #endif ! ***************************************************************** subroutine tmm_Read_TQ( tmm, archivekey_T, archivekey_Q, tday, t1, t2, lli, levi, & sp, T, T_tmi, Q, Q_tmi, status ) use GO , only : goSplitLine use Binas , only : p_global use Phys , only : RealTemperature use tmm_info, only : TMeteoInfo, Init, Done, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey_T character(len=*), intent(in) :: archivekey_Q type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: T(:,:,:) ! K type(TMeteoInfo), intent(out) :: T_tmi real, intent(out) :: Q(:,:,:) ! kg/kg type(TMeteoInfo), intent(out) :: Q_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_TQ' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename ! --- begin ------------------------------- ! split source key in type and name: call goSplitLine( archivekey_T, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) case ( 'standard' ) ! fill field with global mean pressure: sp = p_global ! dummy values: Q = 0.0 T = 0.0 ! set history info call Init( Q_tmi, 'Q', 'kg/kg', status ) call AddHistory( Q_tmi, 'standard', status ) call Init( T_tmi, 'T', 'K', status ) call AddHistory( T_tmi, 'standard', status ) case ( 'tmpp', 'tm5-hdf', 'tm5-nc', 'ecmwf-tmpp', 'ecmwf-tm5', 'msc-data' ) ! humidity: call ReadField( tmm, archivekey_Q, 'Q', 'kg/kg', tday, t1, t2, & lli, 'n', levi, 'n', sp, Q, Q_tmi, status ) IF_NOTOK_RETURN(status=1) ! temperature: call ReadField( tmm, archivekey_T, 'T', 'K', tday, t1, t2, & lli, 'n', levi, 'n', sp, T, T_tmi, status ) IF_NOTOK_RETURN(status=1) case ( 'ncep-cdc', 'ncep-gfs' ) ! humidity: call ReadField( tmm, archivekey_Q, 'Q', 'kg/kg', tday, t1, t2, & lli, 'n', levi, 'n', sp, Q, Q_tmi, status ) IF_NOTOK_RETURN(status=1) ! virtual temperature: call ReadField( tmm, archivekey_T, 'Tv', 'K', tday, t1, t2, & lli, 'n', levi, 'n', sp, T, T_tmi, status ) IF_NOTOK_RETURN(status=1) ! convert: T = RealTemperature( T, Q ) ! info: call AddHistory( T_tmi, 'convert from virtual temperature', status ) case default write (gol,'("unsupported temper source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ! done ! ! ok status = 0 end subroutine tmm_Read_TQ ! ***************************************************************** subroutine tmm_Read_CloudCovers( tmm, archivekey, tday, t1, t2, lli, levi, & sp, cc, cc_tmi, cco, cco_tmi, ccu, ccu_tmi, status ) use GO , only : goSplitLine use Binas , only : p_global use Grid , only : FillLevels use Phys , only : cf_overhead use tmm_info , only : TMeteoInfo, Init, AddHistory use tmm_param, only : GetCombineKeys ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: cc(:,:,:), cco(:,:,:), ccu(:,:,:) ! 0-1 type(TMeteoInfo), intent(out) :: cc_tmi, cco_tmi, ccu_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_CloudCovers' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: i, j, l, lme real, allocatable :: cc_col(:) real, allocatable :: cco_col(:), ccoX(:,:,:) real, allocatable :: ccu_col(:), ccuX(:,:,:) character(len=10) :: hcomb, vcomb ! --- begin ------------------------------- call goLabel(rname) ! check ... if ( any(shape(cco)/=shape(cc)) .or. any(shape(ccu)/=shape(cc)) ) then write (gol,'("output arrays should have same shape:")'); call goErr write (gol,'(" cc : ",3i4)') shape( cc); call goErr write (gol,'(" cco : ",3i4)') shape(cco); call goErr write (gol,'(" ccu : ",3i4)') shape(ccu); call goErr TRACEBACK; status=1; return end if ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! standard values ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'standard' ) ! fill field with global mean pressure: sp = p_global ! no clouds ... cc = 0.0 cco = 0.0 ccu = 0.0 ! set history info call Init( cc_tmi, 'CC', '0-1', status ) call AddHistory( cc_tmi, 'standard', status ) call Init( cco_tmi, 'CCO', '0-1', status ) call AddHistory( cc_tmi, 'standard', status ) call Init( ccu_tmi, 'CCU', '0-1', status ) call AddHistory( cc_tmi, 'standard', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tmpp', 'tm5-hdf', 'tm5-nc', 'prism' ) ! cloud cover call ReadField( tmm, archivekey, 'CC', '1', tday, t1, t2, & lli, 'n', levi, 'n', sp, cc, cc_tmi, status ) IF_NOTOK_RETURN(status=1) ! cloud cover overhead call ReadField( tmm, archivekey, 'CCO', '1', tday, t1, t2, & lli, 'n', levi, 'n', sp, cco, cco_tmi, status ) IF_NOTOK_RETURN(status=1) ! cloud cover underfeet call ReadField( tmm, archivekey, 'CCU', '1', tday, t1, t2, & lli, 'n', levi, 'n', sp, ccu, ccu_tmi, status ) IF_NOTOK_RETURN(status=1) ! set extrema; stored values could be slightly outside [0,1] cc = max( 0.0, min( cc, 1.0 ) ) cco = max( 0.0, min( cco, 1.0 ) ) ccu = max( 0.0, min( ccu, 1.0 ) ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! convert from raw gg fields ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tmpp', 'ecmwf-tm5' ) !DBG call wrtgol( ' tmm read : cc ', t1, ' - ', t2 ); call goPr ! ! read cc: gg, all levels ! call FillBuffer( tmm, archivekey, 'CC', '1', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! 3d horizontal conversion; result is stored in tmm%llX call Transform3Dh( tmm, lli, 'n', levi, 'n', sp, cc_tmi, status ) IF_NOTOK_RETURN(status=1) ! ! compute cloudcover overhead using all levels ! ! unreduced number of levels: lme = tmm%buf_levi%nlev ! storage: allocate( cc_col(lme) ) allocate( cco_col(lme) ) allocate( ccoX(lli%nlon,lli%nlat,lme) ) allocate( ccu_col(lme) ) allocate( ccuX(lli%nlon,lli%nlat,lme) ) ! loop over grid cells do j = 1, lli%nlat do i = 1, lli%nlon ! overhead cloud cover: cc_col = tmm%llX(i,j,:) call cf_overhead ( lme, cc_col, cco_col ) ccoX(i,j,:) = cco_col ! underfeet cloud cover; first reverse layers: do l = 1, lme cc_col(l) = tmm%llX(i,j,lme+1-l) end do call cf_overhead( lme, cc_col, ccu_col ) do l = 1, lme ccuX(i,j,l) = ccu_col(lme+1-l) end do end do end do ! clear deallocate( cc_col ) deallocate( cco_col ) deallocate( ccu_col ) ! info on this operation: call AddHistory( cco_tmi, 'oper==cf_overhead', status ) ! ! 3d vertical conversions ! ! convert from tmm%buf_llX to cc call Transform3Dv( tmm, levi, 'n', sp, cc, cc_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ccoX in buffer, convert from tmm%llX to cco call GetCombineKeys( hcomb, vcomb, 'cco', status ) IF_NOTOK_RETURN(status=1) call FillLevels( levi, 'n', sp, cco, tmm%buf_levi, ccoX, vcomb, status ) IF_NOTOK_RETURN(status=1) call AddHistory( cco_tmi, 'oper==vcomb,'//trim(vcomb), status ) ! store ccuX in buffer, convert from tmm%llX to ccu call GetCombineKeys( hcomb, vcomb, 'ccu', status ) IF_NOTOK_RETURN(status=1) call FillLevels( levi, 'n', sp, ccu, tmm%buf_levi, ccuX, vcomb, status ) IF_NOTOK_RETURN(status=1) call AddHistory( ccu_tmi, 'oper==vcomb,'//trim(vcomb), status ) ! clear deallocate( ccoX ) deallocate( ccuX ) ! set extrema; stored values could be slightly outside [0,1] cc = max( 0.0, min( cc, 1.0 ) ) cco = max( 0.0, min( cco, 1.0 ) ) ccu = max( 0.0, min( ccu, 1.0 ) ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! error ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ok call goLabel() status = 0 end subroutine tmm_Read_CloudCovers ! ***************************************************************** subroutine tmm_Read_Convec( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & gph, gph_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) use GO , only : operator(-), operator(+), rTotal use GO , only : goSplitLine, goVarValue use Binas , only : p_global use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(in) :: omega(:,:,:) ! Pa/s type(TMeteoInfo), intent(in) :: omega_tmi real, intent(in) :: gph(:,:,:) ! m type(TMeteoInfo), intent(in) :: gph_tmi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: entu(:,:,:), entd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: entu_tmi, entd_tmi real, intent(out) :: detu(:,:,:), detd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: detu_tmi, detd_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Convec' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: lout real, allocatable :: ll(:,:,:) character(len=8) :: method ! --- begin ------------------------------- call goLabel(rname) !DBG call wrtgol( ' tmm read convec ', t1, ' - ', t2 ); call goPr ! number of levels in output arrays: lout = size(entu,3) ! check ... if ( (size(entd,3)/=lout) .or. & (size(detu,3)/=lout) .or. (size(detd,3)/=lout) ) then write (gol,'("output arrays should have same number of levels:")'); call goErr write (gol,'(" entu : ",i4)') size(entu,3); call goErr write (gol,'(" entd : ",i4)') size(entd,3); call goErr write (gol,'(" detu : ",i4)') size(detu,3); call goErr write (gol,'(" detd : ",i4)') size(detd,3); call goErr TRACEBACK; status=1; return end if ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! standard values ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'standard') ! fill field with global mean pressure: sp = p_global ! no convection ... entu = 0.0 entd = 0.0 detu = 0.0 detd = 0.0 ! set history info call Init( entu_tmi, 'eu', 'kg/m2/s', status ) call AddHistory( entu_tmi, 'standard', status ) call Init( entd_tmi, 'ed', 'kg/m2/s', status ) call AddHistory( entd_tmi, 'standard', status ) call Init( detu_tmi, 'du', 'kg/m2/s', status ) call AddHistory( detu_tmi, 'standard', status ) call Init( entd_tmi, 'dd', 'kg/m2/s', status ) call AddHistory( detd_tmi, 'standard', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tmpp', 'tm5-hdf', 'tm5-nc' ) ! full level output: allocate( ll(lli%nlon,lli%nlat,levi%nlev) ) ! entrement updraught call ReadField( tmm, archivekey, 'eu', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entu_tmi, status ) IF_NOTOK_RETURN(status=1) ! entu = ll(:,:,1:lout) ! entrement downdraught call ReadField( tmm, archivekey, 'ed', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entd_tmi, status ) IF_NOTOK_RETURN(status=1) ! entd = ll(:,:,1:lout) ! detrement updraught call ReadField( tmm, archivekey, 'du', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detu_tmi, status ) IF_NOTOK_RETURN(status=1) ! detu = ll(:,:,1:lout) ! detrement downdraught call ReadField( tmm, archivekey, 'dd', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detd_tmi, status ) IF_NOTOK_RETURN(status=1) ! detd = ll(:,:,1:lout) ! clear deallocate( ll ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ecmwf convective stuff ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tm5','prism' ) method = 'raw' #ifdef with_prism method = 'ec-ll' #endif ! overwrite method if provided in "source" call goVarValue( sourcename, ';', 'method', '=', method, status ) IF_ERROR_RETURN(status=1) select case ( method ) #ifdef with_tmm_convec_raw case ( 'raw' ) call tmm_Read_Convec_raw( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) IF_NOTOK_RETURN(status=1) #endif #ifdef with_tmm_convec_ec_gg case ( 'ec-gg' ) ! read ec flux/detr, convert to tm entr/detr, average to tm ll call tmm_Read_Convec_EC_gg( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) IF_NOTOK_RETURN(status=1) #endif #ifdef with_tmm_convec_ec case ( 'ec-ll' ) ! read ec flux/detr, aver to tm ll, convert to tm entr/detr ! note: gph instead of omega call tmm_Read_Convec_EC( tmm, archivekey, tday, t1, t2, lli, levi, & gph, gph_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) IF_NOTOK_RETURN(status=1) #endif case default write (gol,'("unsupported convec method : ",a)') trim(method); call goErr TRACEBACK; status=1; return end select #ifdef with_tmm_convec_raw ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! convert from raw fields (sh,gg) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tmpp', 'ncep-cdc', 'ncep-gfs', 'msc-data' ) call tmm_Read_Convec_raw( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) IF_NOTOK_RETURN(status=1) #endif ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! error ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ok call goLabel() status = 0 end subroutine tmm_Read_Convec ! ***************************************************************** subroutine tmm_Read_Diffus( tmm, archivekey, tday, t1, t2, lli, levi, & sp, Kzz, Kzz_tmi, & status ) use GO , only : operator(-), operator(+), rTotal use GO , only : goSplitLine, goVarValue use Binas , only : p_global use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: Kzz(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: Kzz_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Diffus' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: lout real, allocatable :: ll(:,:,:) character(len=8) :: method ! --- begin ------------------------------- call goLabel(rname) !DBG call wrtgol( ' tmm read diffus ', t1, ' - ', t2 ); call goPr ! number of levels in output arrays: lout = size(Kzz,3) ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! standard values ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'standard') ! fill field with global mean pressure: sp = p_global ! no diffusion ... Kzz = 0.0 ! set history info call Init( Kzz_tmi, 'K', 'kg/m2/s', status ) call AddHistory( Kzz_tmi, 'standard', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tmpp', 'tm5-hdf', 'tm5-nc' ) ! half level output: allocate( ll(lli%nlon,lli%nlat,levi%nlev+1) ) ! diffusion flux at half levels: call ReadField( tmm, archivekey, 'K', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'w', sp, ll, Kzz_tmi, status ) IF_NOTOK_RETURN(status=1) ! Kzz = ll(:,:,1:lout) ! clear deallocate( ll ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ecmwf convective stuff ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tm5','prism' ) ! read ec flux/detr, aver to tm ll, convert to tm entr/detr ! note: gph instead of omega call tmm_Read_Diffus_EC( tmm, archivekey, tday, t1, t2, lli, levi, & sp, & Kzz, Kzz_tmi, & status ) IF_NOTOK_RETURN(status=1) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! error ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ok call goLabel() status = 0 end subroutine tmm_Read_Diffus ! ***************************************************************** #ifdef with_tmm_convec_raw subroutine tmm_Read_Convec_raw( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) use Binas , only : grav use GO , only : Get, IncrDate use GO , only : operator(-), operator(+), rTotal, iTotal use GO , only : goSplitLine use Phys , only : mid2bound_uv, mid2bound_w, mid2bound_t use Phys , only : mid2bound_q, mid2bound_z, mid2bound_p use Phys , only : subscal, phlev_ec_rev, geopot use Phys , only : subscal_2d use Phys , only : RealTemperature use Phys , only : GeoPotentialHeightB use Phys , only : ECconv_to_TMconv use Grid , only : HPressure, FPressure, FillLevels use Grid , only : TShGrid, VoD2UV use Grid , only : InterpolMask, Divergence, Tgg2llFracs, FracSum, AreaOper use Grid , only : Init, Done, Set, Interpol use Grid , only : NewInterpol use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(in) :: omega(:,:,:) ! Pa/s type(TMeteoInfo), intent(in) :: omega_tmi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: entu(:,:,:), entd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: entu_tmi, entd_tmi real, intent(out) :: detu(:,:,:), detd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: detu_tmi, detd_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Convec_raw' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: lout real, pointer :: ll(:,:,:) ! timing integer :: hour real :: dhour integer :: dth type(TDate) :: t1s, t2s type(TDate) :: t1m, t2m logical :: skip_second ! loops etc integer :: l, nlev integer :: i, i1, i2, j ! spectral fields type(TshGridInfo) :: shi complex, pointer :: LNSP_sh(:) complex, pointer :: D_sh(:,:), VO_sh(:,:) complex, pointer :: U_sh(:), V_sh(:) ! gaussian grids integer :: ggN logical :: reduced type(TggGridInfo) :: ggi real, pointer :: gg(:) real, pointer :: slhf_gg(:) real, pointer :: oro_gg(:) real, pointer :: sp_gg(:) real, pointer :: m_gg(:) real, pointer :: u_gg(:,:) real, pointer :: v_gg(:,:) real, pointer :: w_gg(:,:) real, pointer :: t_gg(:,:) real, pointer :: q_gg(:,:) real, pointer :: p_gg(:,:) real, pointer :: z_gg(:,:) logical, pointer :: ggmask(:) real, pointer :: qam_gg(:,:), qac_gg(:,:) real, pointer :: entu_gg(:,:) real, pointer :: detu_gg(:,:) real, pointer :: entd_gg(:,:) real, pointer :: detd_gg(:,:) type(Tgg2llFracs) :: gg2ll real, pointer :: llX(:,:,:) ! subgrid stuff real :: dt_sec real :: evap real, pointer :: p_hlev(:) ! extra info type(TMeteoInfo) :: slhf_tmi, sp_tmi, oro_tmi type(TMeteoInfo) :: vo_tmi, div_tmi, u_tmi, v_tmi type(TMeteoInfo) :: w_tmi, t_tmi, q_tmi ! reversed levels type(TLevelInfo) :: leviX real, pointer :: aX(:), bX(:) real, pointer :: tmp_gg(:,:) ! --- begin ------------------------------- ! number of levels in output arrays: lout = size(entu,3) ! check ... if ( (size(entd,3)/=lout) .or. & (size(detu,3)/=lout) .or. (size(detd,3)/=lout) ) then write (gol,'("output arrays should have same number of levels:")'); call goErr write (gol,'(" entu : ",i4)') size(entu,3); call goErr write (gol,'(" entd : ",i4)') size(entd,3); call goErr write (gol,'(" detu : ",i4)') size(detu,3); call goErr write (gol,'(" detd : ",i4)') size(detd,3); call goErr TRACEBACK; status=1; return end if ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! length of time interval: dth = iTotal( t2 - t1, 'hour' ) dt_sec = dth * 3600.0 ! 3 hourly or 6 hourly ? select case ( dth ) ! ! ~~ 6 hourly ! case ( 6 ) ! only around hours 00/06/12/ etc yet ... call Get( t1, hour=hour ) if ( modulo(hour,6) /= 3 ) then write (gol,'("6 hourly convection only for intervals [21,03], [03,09], ...")'); call goErr TRACEBACK; status=1; return end if ! times for model level fields is mid of requested interval: t1m = t1 + IncrDate(sec=nint(dt_sec/2)) t2m = t1m ! use 6 hour interval around requested (instant!) time to read slhf; ! for 3 hourly request, this means double reading ! (no problem since slhf is time average anyway) ! ! slhf (gg) ! ! first interval: [t1-dt/2,t1] t1s = t1 - IncrDate(sec=int(dt_sec/2)) t2s = t1 !call wrtgol( ' tmm read : slhf ', t1s, ' - ', t2s ); call goPr ! read first slhf in buffer (W/m2 time aver): call FillBuffer( tmm, archivekey, 'slhf', 'W/m2', tday, t1s, t2s, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg' ) IF_NOTOK_RETURN(status=1) ! extract grid size ggN = tmm%buf_ggi%N reduced = tmm%buf_ggi%reduced ! setup gg defintion from buffer: call Init( ggi, ggN, reduced, status ) IF_NOTOK_RETURN(status=1) ! allocate storage: allocate( oro_gg(ggi%np) ) allocate( slhf_gg(ggi%np) ) allocate( sp_gg(ggi%np) ) ! store first half of slhf (accumulated over 3 hr): slhf_gg = tmm%buf_gg(:,1) * dt_sec/2 ! W/m2 s slhf_tmi = tmm%buf_tmi ! second interval: [t1,t1+dt/2] t1s = t1 t2s = t1 + IncrDate(sec=int(dt_sec/2)) !call wrtgol( ' tmm read : slhf ', t1s, ' - ', t2s ); call goPr ! read first slhf in buffer (W/m2 time aver): call FillBuffer( tmm, archivekey, 'slhf', 'W/m2', tday, t1s, t2s, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np ) IF_NOTOK_RETURN(status=1) ! add second half of slhf (accumulated over 3 hr): slhf_gg = slhf_gg + tmm%buf_gg(:,1) * dt_sec/2 ! W/m2 s ! copy info call AddHistory( slhf_tmi, tmm%buf_tmi, status ) ! ! ~~ 3 hourly ! #ifdef with_tmm_ecearth case ( 0 ) ! assume this is the 6 hourly meteo ... ! only for times around 00, 06, ... call Get( t1, hour=hour ) if ( modulo(hour,6) /= 0 ) then write (gol,'("6 hourly convection only for intervals [00,06], [06,12], ...")'); call goErr TRACEBACK; status=1; return end if ! times for model level fields is begin of requested interval t1m = t1 t2m = t1m ! use 12 hour interval around requested (instant!) time to read slhf; ! ! slhf (gg), accumulated over [-6,6] hours ! ! slhf will be valid for [-6,6] dt_sec = 12*3600.0 ! second field should be read skip_second = .false. #else case ( 0, 3 ) ! convective stuff for instant time; ! assume this is the 3 hourly meteo ... ! only for times around 00, 03, ... call Get( t1, hour=hour ) if ( modulo(hour,3) /= 0 ) then write (gol,'("3 hourly convection only for intervals [00,03], [03,06], ...")'); call goErr TRACEBACK; status=1; return end if ! times for model level fields is begin of requested interval ! (just a choice) t1m = t1 t2m = t1m ! use 6 hour interval around requested (instant!) time to read slhf; ! for 3 hourly request, this means double reading ! (no problem since slhf is time average anyway) ! ! slhf (gg), accumulated over [-3,3] hours ! ! slhf will be valid for [-3,3] dt_sec = 6*3600.0 ! forecast for 72 hour or later ? then slhf for [-6,6] if ( rTotal(t1-tday,'hour') >= 12.0+72.0 ) dt_sec = 12*3600.0 ! forecast for fcday 10 is only for [00:00,12:00] where 12:00 is 12+240, ! thus no second interval; instead copy from first skip_second = rTotal(t1-tday,'hour') >= 12.0+240.0 ! >>> adhoc: always skip second interval, ! because of problems with setting tday; ! for requested time, this could be previous day ! while slhf is taken from two different days ! (in fact routines should be called with two tdays rather than one) !skip_second = .true. !write (gol,'("WARNING - skipped second interval for slhf")'); call goPr #endif ! first interval: [t1-dt,t1] t1s = t1 - IncrDate(sec=int(dt_sec/2)) t2s = t1 !call wrtgol( ' tmm read : slhf ', t1s, ' - ', t2s ); call goPr ! read first slhf in buffer (W/m2 time aver): call FillBuffer( tmm, archivekey, 'slhf', 'W/m2', tday, t1s, t2s, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg' ) IF_NOTOK_RETURN(status=1) ! extract grid size ggN = tmm%buf_ggi%N reduced = tmm%buf_ggi%reduced ! setup gg defintion from buffer: call Init( ggi, ggN, reduced, status ) IF_NOTOK_RETURN(status=1) ! allocate storage: allocate( oro_gg(ggi%np) ) allocate( slhf_gg(ggi%np) ) allocate( sp_gg(ggi%np) ) ! store first half of slhf (accumulated over 3/6 hr): slhf_gg = tmm%buf_gg(:,1) * dt_sec/2 ! W/m2 s slhf_tmi = tmm%buf_tmi ! second interval: [t1,t1+dt] t1s = t1 t2s = t1 + IncrDate(sec=int(dt_sec/2)) ! skip seond interval ? if ( skip_second ) then !DBG call wrtgol( 'tmm copy : slhf ', t1s, ' - ', t2s ); call goPr ! use for the second field the first: slhf_gg = slhf_gg * 2 else !call wrtgol( ' tmm read : slhf ', t1s, ' - ', t2s ); call goPr ! read second slhf in buffer (W/m2 time aver): call FillBuffer( tmm, archivekey, 'slhf', 'W/m2', tday, t1s, t2s, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np ) IF_NOTOK_RETURN(status=1) ! add second half of slhf (accumulated over 3 hr): slhf_gg = slhf_gg + tmm%buf_gg(:,1) * dt_sec/2 ! W/m2 s ! copy info call AddHistory( slhf_tmi, tmm%buf_tmi, status ) end if ! ! ~~ error ! case default write (gol,'("unsupported lenght of time interval : ",i4)') dth; call goErr TRACEBACK; status=1; return end select ! ! orography (gg) ! !write (gol,'(" tmm read : oro")'); call goPr ! read orography in buffer: call FillBuffer( tmm, archivekey, 'oro', 'm m/s2', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! interpol from sh ? select case ( tmm%buf_gridtype ) case ( 'gg' ) ! copy from buffer: oro_gg = tmm%buf_gg(:,1) ! m*[g] case ( 'sh' ) ! interpol from sh to gg : call Interpol( tmm%buf_shi, tmm%buf_sh(:,1), ggi, oro_gg, status ) IF_NOTOK_RETURN(status=1) case default write (gol,'("unsupported grid type `",a,"` for raw oro")') tmm%buf_gridtype TRACEBACK; status=1; return end select ! store: oro_tmi = tmm%buf_tmi ! ! read raw 3d fields ! ! ~~~ !call wrtgol( ' tmm read : u,v ', t1m, ' - ', t2m ); call goPr ! read VO call FillBuffer( tmm, archivekey, 'VO', '1/s', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh' ) IF_NOTOK_RETURN(status=1) ! extract other grid sizes call Init( shi, tmm%buf_shi, status ) IF_NOTOK_RETURN(status=1) nlev = tmm%buf_levi%nlev ! allocate 3d storage: allocate( u_gg(ggi%np,0:nlev) ) allocate( v_gg(ggi%np,0:nlev) ) allocate( w_gg(ggi%np,0:nlev) ) allocate( t_gg(ggi%np,0:nlev) ) allocate( q_gg(ggi%np,0:nlev) ) allocate( p_gg(ggi%np,0:nlev) ) allocate( z_gg(ggi%np,0:nlev) ) ! copy 3d spectral field from buffer: allocate( VO_sh(shi%np,nlev) ) VO_sh = tmm%buf_sh ! copy info vo_tmi = tmm%buf_tmi ! read D call FillBuffer( tmm, archivekey, 'D', '1/s', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d spectral field from buffer: allocate( D_sh(shi%np,nlev) ) D_sh = tmm%buf_sh ! copy info div_tmi = tmm%buf_tmi ! convert from sh VO/D to gg U/V ! setup storage: allocate( U_sh(shi%np) ) allocate( V_sh(shi%np) ) ! loop over levels: !xOMP PARALLEL & !xOMP default ( none ) & !xOMP shared ( nlev, shi, VO_sh, D_sh ) & !xOMP shared ( ggi, u_gg, v_gg ) & !xOMP private ( l, U_sh, V_sh, status ) !xOMP DO do l = 1, nlev ! convert to U=u*cos(lat) and V=v*cos(lat) : call vod2uv( shi, VO_sh(:,l), D_sh(:,l), shi, U_sh, V_sh ) ! convert to Gaussian grid: call NewInterpol( shi, U_sh, ggi, u_gg(:,l), status ) !IF_NOTOK_RETURN(status=1) call NewInterpol( shi, V_sh, ggi, v_gg(:,l), status ) !IF_NOTOK_RETURN(status=1) end do !xOMP END DO !xOMP END PARALLEL ! clear call Done( shi ) deallocate( U_sh ) deallocate( V_sh ) ! history ... call Init( u_tmi, 'u', 'm/s', status, (/vo_tmi,div_tmi/) ) call Init( v_tmi, 'u', 'm/s', status, (/vo_tmi,div_tmi/) ) call AddHistory( u_tmi, 'VoD to UV;;interpol to gg', status ) call AddHistory( v_tmi, 'VoD to UV;;interpol to gg', status ) ! clear deallocate( VO_sh ) deallocate( D_sh ) ! remove cos(lat) factor: do j = 1, ggi%nlat i1 = ggi%i1(j) i2 = ggi%im(j) u_gg(i1:i2,1:nlev) = u_gg(i1:i2,1:nlev) / cos(ggi%lat(j)) v_gg(i1:i2,1:nlev) = v_gg(i1:i2,1:nlev) / cos(ggi%lat(j)) end do call AddHistory( u_tmi, 'divide by cos(lat)', status ) call AddHistory( v_tmi, 'divide by cos(lat)', status ) ! ~~~ select case ( sourcetype ) case ( 'ecmwf-tmpp', 'msc-data' ) !call wrtgol( ' tmm read : Q ', t1m, ' - ', t2m ); call goPr ! read humidity call FillBuffer( tmm, archivekey, 'Q', 'kg/kg', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d gg field from buffer; copy info: Q_gg(:,1:nlev) = tmm%buf_gg ! kg h2o / kg air ! info ... call Init( Q_tmi, tmm%buf_tmi, status ) ! copy surface pressure: sp_gg = tmm%buf_sp_gg ! Pa call Init( sp_tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) !call wrtgol( ' tmm read : T ', t1m, ' - ', t2m ); call goPr ! read T call FillBuffer( tmm, archivekey, 'T', 'K', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T, nlev=nlev ) IF_NOTOK_RETURN(status=1) do l = 1, nlev call Interpol( tmm%buf_shi, tmm%buf_sh(:,l), ggi, T_gg(:,l), status ) IF_NOTOK_RETURN(status=1) end do ! info ... call Init( T_tmi, tmm%buf_tmi, status ) call AddHistory( T_tmi, 'interpol to gg', status ) case ( 'ecmwf-tm5' ) !call wrtgol( ' tmm read : Q ', t1m, ' - ', t2m ); call goPr ! read humidity call FillBuffer( tmm, archivekey, 'Q', 'kg/kg', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d gg field from buffer:: Q_gg(:,1:nlev) = tmm%buf_gg ! kg h2o / kg air ! info ... call Init( Q_tmi, tmm%buf_tmi, status ) ! copy surface pressure: sp_gg = tmm%buf_sp_gg ! Pa call Init( sp_tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) !call wrtgol( ' tmm read : T ', t1m, ' - ', t2m ); call goPr ! read T call FillBuffer( tmm, archivekey, 'T', 'K', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d gg field from buffer:: T_gg(:,1:nlev) = tmm%buf_gg ! K ! info ... call Init( T_tmi, tmm%buf_tmi, status ) case ( 'ncep-cdc', 'ncep-gfs' ) !call wrtgol( ' tmm read : Q ', t1m, ' - ', t2m ); call goPr ! read humidity call FillBuffer( tmm, archivekey, 'Q', 'kg/kg', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', np=ggi%np, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! convert from sh to gg do l = 1, nlev call Interpol( tmm%buf_shi, tmm%buf_sh(:,l), ggi, Q_gg(:,l), status ) ! kg h2o / kg air IF_NOTOK_RETURN(status=1) end do ! prevent negatives ... Q_gg = max( 0.0, Q_gg ) ! info ... call Init( Q_tmi, tmm%buf_tmi, status ) call AddHistory( Q_tmi, 'interpol to gg', status ) call AddHistory( Q_tmi, 'truncate', status ) ! interpolate surface pressure: call Interpol( tmm%buf_shi, tmm%buf_lnsp_sh, ggi, sp_gg, status ) ! ln(Pa) IF_NOTOK_RETURN(status=1) sp_gg = exp( sp_gg ) ! Pa ! info ... call Init( sp_tmi, 'sp', 'Pa', status, (/tmm%buf_tmi/) ) !call wrtgol( ' tmm read : T ', t1m, ' - ', t2m ); call goPr ! read virtual temperature call FillBuffer( tmm, archivekey, 'Tv', 'K', tday, t1m, t2m, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='sh', shT=shi%T, nlev=nlev ) IF_NOTOK_RETURN(status=1) ! convert from sh to gg do l = 1, nlev call Interpol( tmm%buf_shi, tmm%buf_sh(:,l), ggi, T_gg(:,l), status ) IF_NOTOK_RETURN(status=1) end do ! convert from virtual to normal temperature: T_gg = RealTemperature( T_gg, Q_gg ) ! info ... call Init( T_tmi, tmm%buf_tmi, status ) call AddHistory( T_tmi, 'interpol to gg', status ) case default write (gol,'("unsupported source type `",a,"` for raw surface fields")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ~~~ write (gol,'(" interpol : W")'); call goPr ! omega on gg and TM levels: allocate( tmp_gg(ggi%np,levi%nlev+1) ) ! convert from ll to gg: do l = 1, levi%nlev+1 call Interpol( lli, omega(:,:,l), ggi, tmp_gg(:,l), status ) IF_NOTOK_RETURN(status=1) end do ! convert to parent levels (buffer) from TM levels: call FillLevels( tmm%buf_levi, 'w', sp_gg, w_gg, & levi, tmp_gg, 'bottom', status ) IF_NOTOK_RETURN(status=1) ! info ... call Init( w_tmi, omega_tmi, status ) call AddHistory( w_tmi, 'interpol to gg and raw levels', status ) ! clear deallocate( tmp_gg ) ! ! ensure that layers are in ecmwf order (top->down) ! select case ( sourcetype ) case ( 'ecmwf-tmpp', 'ecmwf-tm5', 'msc-data' ) ! copy level info from buffer: call Init( leviX, tmm%buf_levi%nlev, tmm%buf_levi%a, tmm%buf_levi%b, status ) IF_NOTOK_RETURN(status=1) case ( 'ncep-cdc', 'ncep-gfs' ) ! revert level info from buffer: allocate( aX(0:tmm%buf_levi%nlev) ) allocate( bX(0:tmm%buf_levi%nlev) ) do l = 0, tmm%buf_levi%nlev aX(l) = tmm%buf_levi%a(tmm%buf_levi%nlev-l) bX(l) = tmm%buf_levi%b(tmm%buf_levi%nlev-l) end do call Init( leviX, tmm%buf_levi%nlev, aX, bX, status ) IF_NOTOK_RETURN(status=1) deallocate( aX ) deallocate( bX ) ! revert fields allocate( tmp_gg(ggi%np,0:nlev) ) tmp_gg = u_gg do l = 1, nlev u_gg(:,l) = tmp_gg(:,nlev+1-l) end do tmp_gg = v_gg do l = 1, nlev v_gg(:,l) = tmp_gg(:,nlev+1-l) end do tmp_gg = w_gg do l = 1, nlev w_gg(:,l) = tmp_gg(:,nlev+1-l) end do tmp_gg = t_gg do l = 1, nlev t_gg(:,l) = tmp_gg(:,nlev+1-l) end do tmp_gg = q_gg do l = 1, nlev q_gg(:,l) = tmp_gg(:,nlev+1-l) end do deallocate( tmp_gg ) end select ! ! updraughts/downdraughts ! write (gol,'(" tmm updr/downdr")'); call goPr ! ! WARNING: the TMPP2/tmpp_conv-gg on bsgi59 is probably not bugfree: ! o oro not read ! o pw = w/g [kg/m2/s], but tmpp_conv_tiedtke expects Pa/s ? ! o factor to convert from lshf to evaporation ! here : 1 m/s = 2.45e9 W/m2 (at 20 degrees C) ! binas : lvap = 2.5 e6 J/kg ! binas : Lc = 2.26e6 J/kg (at 0 deg C) ! ! There are however good reasons to step over to the conv-gg code ! by Dirk Olivie ! o no unnecessary interpolations to half levels ! o removed stuff that was based on the ec 19 levels ! ! extra fields: allocate( ggmask(ggi%np) ) allocate( qam_gg(ggi%np,0:nlev) ) allocate( qac_gg(ggi%np,0:nlev) ) allocate( entu_gg(ggi%np,nlev) ) allocate( detu_gg(ggi%np,nlev) ) allocate( entd_gg(ggi%np,nlev) ) allocate( detd_gg(ggi%np,nlev) ) ! Set mask for averaging over ll; ! one extra row to compute derivatives. ! This routine changes ggi: flag for each row to be processed or not call InterpolMask( ggi, ggmask, lli, 1 ) ! calculate geopotential (z) allocate( p_hlev(0:nlev) ) do i = 1, ggi%np ! skip if not required for ll grid: if ( .not. ggmask(i) ) cycle ! compute pressure at half levels (surf -> top): call phlev_ec_rev( nlev, leviX%a, leviX%b, sp_gg(i), p_hlev ) ! compute z for single column: call GeoPot( nlev, oro_gg(i), T_gg(i,1:nlev), Q_gg(i,1:nlev), & p_hlev, z_gg(i,1:nlev) ) end do deallocate( p_hlev ) ! interpolate variables from the center of parent-model layers to the ! boundaries of parent-model layers and save result in same memory location .. call mid2bound_uv( nlev, ggi%np, u_gg, v_gg, sp_gg, ggmask, leviX%a, leviX%b ) call mid2bound_t ( nlev, ggi%np, t_gg, sp_gg, ggmask, leviX%a, leviX%b ) call mid2bound_q ( nlev, ggi%np, q_gg, sp_gg, ggmask, leviX%a, leviX%b, t_gg ) call mid2bound_z ( nlev, ggi%np, z_gg, sp_gg, ggmask, leviX%a, leviX%b, oro_gg ) ! already on half levels, since interpolated from omega ... !!call mid2bound_w ( nlev, ggi%np, w_gg, sp_gg, ggmask, leviX%a, leviX%b ) ! NOTE: p is not filled on input, but filled with pressures on output call mid2bound_p ( nlev, ggi%np, p_gg, sp_gg, ggmask, leviX%a, leviX%b ) ! divergence fields do l = 0, nlev call Divergence( ggi, q_gg(:,l)*u_gg(:,l), q_gg(:,l)*v_gg(:,l), qac_gg(:,l) ) call Divergence( ggi, u_gg(:,l), v_gg(:,l), qam_gg(:,l) ) end do ! Convert from SLHF (W/m2*interval) to EVAP (m/s) ! 1 m/s = 2.45e9 W/m2 (at 20 degrees C). ! Don't forget to change sign from latent heatflux to evaporation!!! ! evap = - slhf_gg(i) / dt_sec / 2.45e9 ! m/s ! (apply in argument) ! work routine: call subscal_2d( ggi%np, nlev, leviX%a, leviX%b, & z_gg, p_gg, w_gg, t_gg, & q_gg, qac_gg, qam_gg, -1.0e3*slhf_gg/dt_sec/2.45e9, & entu_gg, detu_gg, entd_gg, detd_gg ) ! clear deallocate( ggmask ) deallocate( slhf_gg ) deallocate( oro_gg ) deallocate( u_gg ) deallocate( v_gg ) deallocate( w_gg ) deallocate( t_gg ) deallocate( q_gg ) deallocate( qam_gg ) deallocate( qac_gg ) deallocate( p_gg ) deallocate( z_gg ) ! history call Init( entu_tmi, 'entu', 'kg/m2/s', status, & (/sp_tmi,slhf_tmi,oro_tmi,T_tmi,Q_tmi,u_tmi,v_tmi,w_tmi/) ) call Init( entd_tmi, 'entd', 'kg/m2/s', status, & (/sp_tmi,slhf_tmi,oro_tmi,T_tmi,Q_tmi,u_tmi,v_tmi,w_tmi/) ) call Init( detu_tmi, 'detu', 'kg/m2/s', status, & (/sp_tmi,slhf_tmi,oro_tmi,T_tmi,Q_tmi,u_tmi,v_tmi,w_tmi/) ) call Init( detd_tmi, 'detd', 'kg/m2/s', status, & (/sp_tmi,slhf_tmi,oro_tmi,T_tmi,Q_tmi,u_tmi,v_tmi,w_tmi/) ) ! ! convert from gg to ll ! ! determine fraction call Init( gg2ll, ggi, lli, status ) IF_NOTOK_RETURN(status=1) ! take fractions of overlapping cells call FracSum( gg2ll, sp_gg, sp, status, 'area-aver' ) IF_NOTOK_RETURN(status=1) ! clear deallocate( sp_gg ) ! full level output: allocate( llX(lli%nlon,lli%nlat,nlev ) ) allocate( ll (lli%nlon,lli%nlat,levi%nlev) ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, entu_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: entu = ll(:,:,1:lout) ! info .. call AddHistory( entu_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, entd_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: entd = ll(:,:,1:lout) ! info .. call AddHistory( entd_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, detu_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: detu = ll(:,:,1:lout) ! info .. call AddHistory( detu_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, detd_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: detd = ll(:,:,1:lout) ! info .. call AddHistory( detd_tmi, 'gg to ll, area aver;;sum levels', status ) ! clear call Done( gg2ll ) deallocate( entu_gg ) deallocate( entd_gg ) deallocate( detu_gg ) deallocate( detd_gg ) deallocate( llX ) deallocate( ll ) call Done( leviX, status ) IF_NOTOK_RETURN(status=1) ! ! done ! call Done( ggi, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine tmm_Read_Convec_raw #endif ! ***************************************************************** #ifdef with_tmm_convec_ec_gg subroutine tmm_Read_Convec_EC_gg( tmm, archivekey, tday, t1, t2, lli, levi, & omega, omega_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) use Binas , only : grav use GO , only : Get, IncrDate use GO , only : operator(-), operator(+), rTotal use GO , only : goSplitLine use Phys , only : mid2bound_uv, mid2bound_w, mid2bound_t use Phys , only : mid2bound_q, mid2bound_z, mid2bound_p use Phys , only : subscal, phlev_ec_rev, geopot use Phys , only : RealTemperature use Phys , only : GeoPotentialHeightB use Phys , only : ECconv_to_TMconv use Grid , only : HPressure, FPressure, FillLevels use Grid , only : TShGrid, VoD2UV use Grid , only : InterpolMask, Divergence, Tgg2llFracs, FracSum, AreaOper use Grid , only : Init, Done, Set, Interpol use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(in) :: omega(:,:,:) ! Pa/s type(TMeteoInfo), intent(in) :: omega_tmi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: entu(:,:,:), entd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: entu_tmi, entd_tmi real, intent(out) :: detu(:,:,:), detd(:,:,:) ! kg/m2/s type(TMeteoInfo), intent(out) :: detu_tmi, detd_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Convec_EC_gg' ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: lout real, allocatable :: ll(:,:,:) ! timing integer :: hour real :: dhour type(TDate) :: t1s, t2s ! loops etc integer :: l, nlev integer :: i, i1, i2, j ! gaussian grids integer :: ggN logical :: reduced type(TggGridInfo) :: ggi real, allocatable :: gg(:) real, allocatable :: oro_gg(:) real, allocatable :: sp_gg(:) real, allocatable :: m_gg(:) real, allocatable :: t_gg(:,:) real, allocatable :: q_gg(:,:) type(Tgg2llFracs) :: gg2ll real, allocatable :: llX(:,:,:) type(TLevelInfo) :: leviX real, allocatable :: p_hlev(:) ! e4 convection fields real, allocatable :: udmf_gg(:,:) real, allocatable :: ddmf_gg(:,:) real, allocatable :: uddr_gg(:,:) real, allocatable :: dddr_gg(:,:) type(TMeteoInfo) :: udmf_tmi, ddmf_tmi, uddr_tmi, dddr_tmi real, allocatable :: ph_ec(:) real, allocatable :: zh_ec(:) real, pointer :: entu_gg(:,:) real, pointer :: detu_gg(:,:) real, pointer :: entd_gg(:,:) real, pointer :: detd_gg(:,:) ! --- begin ------------------------------- ! number of levels in output arrays: lout = size(entu,3) ! check ... if ( (size(entd,3)/=lout) .or. & (size(detu,3)/=lout) .or. (size(detd,3)/=lout) ) then write (gol,'("output arrays should have same number of levels:")'); call goErr write (gol,'(" entu : ",i4)') size(entu,3); call goErr write (gol,'(" entd : ",i4)') size(entd,3); call goErr write (gol,'(" detu : ",i4)') size(detu,3); call goErr write (gol,'(" detd : ",i4)') size(detd,3); call goErr TRACEBACK; status=1; return end if ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! ! Original parameters archived in MARS: ! ! Parameter Surfaces Code Units Units ! 1) 2) 3) ! ! updraught mass flux half levels UDMF 101 kg/m2 kg/m2/s ! downdraught mass flux half levels DDMF 102 kg/m2 kg/m2/s ! updraught detrainment rate full levels UDDR 103 kg/m3 kg/m3/s ! downdraught detrainment rate full levels DDDR 104 kg/m3 kg/m3/s ! ! 1) GRIB code table 2 version 128 ! 2) original units, accumulated ! 3) time averaged after reading ! ! only hours 00/03/06/etc yet ... call Get( t1, hour=hour ) dhour = rTotal( t2 - t1, 'hour' ) if ( (modulo(hour,3) /= 0) .or. (dhour /= 3.0) ) then write (gol,'("convection only for 3hr intervals [0,3] etc.")'); call goErr call wrtgol( ' requested : ', t1, ' - ', t2 ); call goErr write (gol,'(" dhour : ",f8.4)') dhour; call goErr TRACEBACK; status=1; return end if ! ! read gg fields ! !call wrtgol( ' tmm read : UDMF ', t1, ' - ', t2 ); call goPr ! read updraught mass flux call FillBuffer( tmm, archivekey, 'UDMF', 'kg/m2/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg' ) IF_NOTOK_RETURN(status=1) ! extract grid sizes ggN = tmm%buf_ggi%N reduced = tmm%buf_ggi%reduced ! copy level info from buffer: call Init( leviX, tmm%buf_levi%nlev, tmm%buf_levi%a, tmm%buf_levi%b, status ) IF_NOTOK_RETURN(status=1) ! setup gg defintion from buffer: call Init( ggi, ggN, reduced, status ) IF_NOTOK_RETURN(status=1) ! allocate storage: allocate( udmf_gg(ggi%np,0:leviX%nlev) ) allocate( ddmf_gg(ggi%np,0:leviX%nlev) ) allocate( uddr_gg(ggi%np,1:leviX%nlev) ) allocate( dddr_gg(ggi%np,1:leviX%nlev) ) allocate( sp_gg(ggi%np) ) allocate( oro_gg(ggi%np) ) allocate( T_gg(ggi%np,1:leviX%nlev) ) allocate( Q_gg(ggi%np,1:leviX%nlev) ) ! copy 3d field, levels top down, surface implicit zero, copy info: udmf_gg(:,0:nlev-1) = tmm%buf_gg ! kg/m2/s udmf_gg(:, nlev) = 0.0 call AddHistory( udmf_tmi, tmm%buf_tmi, status ) !call wrtgol( ' tmm read : DDMF ', t1, ' - ', t2 ); call goPr ! downdraught mass flux call FillBuffer( tmm, archivekey, 'DDMF', 'kg/m2/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d field, levels top down, surface implicit zero, copy info: ddmf_gg(:,0:nlev-1) = tmm%buf_gg ! kg/m2/s ddmf_gg(:, nlev) = 0.0 call AddHistory( ddmf_tmi, tmm%buf_tmi, status ) !call wrtgol( ' tmm read : UDDR ', t1, ' - ', t2 ); call goPr ! updraught detrainment rate call FillBuffer( tmm, archivekey, 'UDDR', 'kg/m3/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d field, levels top down, copy info: uddr_gg = tmm%buf_gg ! kg/m3/s call AddHistory( uddr_tmi, tmm%buf_tmi, status ) !call wrtgol( ' tmm read : DDDR ', t1, ' - ', t2 ); call goPr ! downdraught detrainment rate call FillBuffer( tmm, archivekey, 'DDDR', 'kg/m3/s', tday, t1, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d field, levels top down, copy info: dddr_gg = tmm%buf_gg ! kg/m3/s call AddHistory( dddr_tmi, tmm%buf_tmi, status ) ! temperature at begin of interval !call wrtgol( ' tmm read : T ', t1, ' - ', t1 ); call goPr call FillBuffer( tmm, archivekey, 'T', 'K', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d field: T_gg = tmm%buf_gg ! K ! temperature at end of interval !call wrtgol( ' tmm read : T ', t2, ' - ', t2 ); call goPr call FillBuffer( tmm, archivekey, 'T', 'K', tday, t2, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! add, average: T_gg = ( T_gg + tmm%buf_gg )/2.0 ! K ! humidity at begin of interval !call wrtgol( ' tmm read : Q ', t1, ' - ', t1 ); call goPr call FillBuffer( tmm, archivekey, 'Q', 'kg/kg', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! copy 3d field: Q_gg = tmm%buf_gg ! kg/kg ! humidity at end of interval !call wrtgol( ' tmm read : Q ', t2, ' - ', t2 ); call goPr call FillBuffer( tmm, archivekey, 'Q', 'kg/kg', tday, t2, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np, nlev=leviX%nlev ) IF_NOTOK_RETURN(status=1) ! add, average: Q_gg = ( Q_gg + tmm%buf_gg )/2.0 ! kg/kg ! surface pressure at begin of interval !call wrtgol( ' tmm read : sp ', t1, ' - ', t1 ); call goPr call FillBuffer( tmm, archivekey, 'sp', 'Pa', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np ) IF_NOTOK_RETURN(status=1) ! copy 2d field: sp_gg = tmm%buf_gg(:,1) ! Pa ! surface pressure at end of interval !call wrtgol( ' tmm read : sp ', t2, ' - ', t2 ); call goPr call FillBuffer( tmm, archivekey, 'sp', 'Pa', tday, t2, t2, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np ) IF_NOTOK_RETURN(status=1) ! add, average: sp_gg = ( sp_gg + tmm%buf_gg(:,1) )/2.0 ! Pa ! read orography in buffer: !write (gol,'(" tmm read : oro")'); call goPr call FillBuffer( tmm, archivekey, 'oro', 'm m/s2', tday, t1, t1, 'n', 'n', status ) IF_NOTOK_RETURN(status=1) ! check ... call CheckBuffer( tmm, status, gridtype='gg', np=ggi%np ) IF_NOTOK_RETURN(status=1) ! copy from buffer: oro_gg = tmm%buf_gg(:,1) ! m*[g] ! ! convert ! allocate( ph_ec(0:leviX%nlev) ) allocate( zh_ec(0:leviX%nlev) ) allocate( entu_gg(ggi%np,leviX%nlev) ) allocate( detu_gg(ggi%np,leviX%nlev) ) allocate( entd_gg(ggi%np,leviX%nlev) ) allocate( detd_gg(ggi%np,leviX%nlev) ) ! loop over gg cells: do i = 1, ggi%np ! half level pressure: call HPressure( leviX, sp_gg(i), ph_ec, status ) IF_NOTOK_RETURN(status=1) ! gph at half levels: call GeoPotentialHeightB( nlev, ph_ec, T_gg(i,:), Q_gg(i,:), oro_gg(i)/grav, zh_ec ) ! convert from fluxes to rates: call ECconv_to_TMconv( leviX%nlev, zh_ec, & udmf_gg(i,:), uddr_gg(i,:), ddmf_gg(i,:), dddr_gg(i,:), & entu_gg(i,:), detu_gg(i,:), entd_gg(i,:), detd_gg(i,:), & status ) IF_NOTOK_RETURN(status=1) end do deallocate( ph_ec ) deallocate( zh_ec ) ! clear deallocate( udmf_gg ) deallocate( ddmf_gg ) deallocate( uddr_gg ) deallocate( dddr_gg ) deallocate( oro_gg ) deallocate( T_gg ) deallocate( Q_gg ) ! history call Init( entu_tmi, 'entu', 'kg/m2/s', status, (/udmf_tmi,ddmf_tmi,uddr_tmi,dddr_tmi/) ) call Init( entd_tmi, 'entd', 'kg/m2/s', status, (/udmf_tmi,ddmf_tmi,uddr_tmi,dddr_tmi/) ) call Init( detu_tmi, 'detu', 'kg/m2/s', status, (/udmf_tmi,ddmf_tmi,uddr_tmi,dddr_tmi/) ) call Init( detd_tmi, 'detd', 'kg/m2/s', status, (/udmf_tmi,ddmf_tmi,uddr_tmi,dddr_tmi/) ) ! ! convert from gg to ll ! ! determine fraction call Init( gg2ll, ggi, lli, status ) IF_NOTOK_RETURN(status=1) ! take fractions of overlapping cells call FracSum( gg2ll, sp_gg, sp, status, 'area-aver' ) IF_NOTOK_RETURN(status=1) ! full level output: allocate( llX(lli%nlon,lli%nlat,nlev ) ) allocate( ll (lli%nlon,lli%nlat,levi%nlev) ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, entu_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: entu = ll(:,:,1:lout) ! info .. call AddHistory( entu_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, entd_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: entd = ll(:,:,1:lout) ! info .. call AddHistory( entd_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, detu_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: detu = ll(:,:,1:lout) ! info .. call AddHistory( detu_tmi, 'gg to ll, area aver;;sum levels', status ) ! take fractions of overlapping cells do l = 1, nlev call FracSum( gg2ll, detd_gg(:,l), llX(:,:,l), status, 'area-aver' ) IF_NOTOK_RETURN(status=1) end do ! integrated variables might become slightly negative ... llX = max( 0.0, llX ) ! combine levels from llX to ll: call FillLevels( levi, 'n', sp, ll, leviX, llX, 'sum', status ) IF_NOTOK_RETURN(status=1) ! truncate to number of output levels: detd = ll(:,:,1:lout) ! info .. call AddHistory( detd_tmi, 'gg to ll, area aver;;sum levels', status ) ! clear call Done( gg2ll ) deallocate( entu_gg ) deallocate( entd_gg ) deallocate( detu_gg ) deallocate( detd_gg ) deallocate( llX ) deallocate( ll ) ! clear deallocate( sp_gg ) ! ! done ! call Done( ggi, status ) IF_NOTOK_RETURN(status=1) call Done( leviX, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine tmm_Read_Convec_EC_gg #endif ! ***************************************************************** #ifdef with_tmm_convec_ec ! ! Original parameters archived in MARS: ! ! Parameter Surfaces Code Units Units ! 1) 2) 3) ! ! updraught mass flux half levels UDMF 101 kg/m2 kg/m2/s ! downdraught mass flux half levels DDMF 102 kg/m2 kg/m2/s ! updraught detrainment rate full levels UDDR 103 kg/m3 kg/m3/s ! downdraught detrainment rate full levels DDDR 104 kg/m3 kg/m3/s ! ! 1) GRIB code table 2 version 128 ! 2) original units, accumulated ! 3) time averaged after reading ! subroutine tmm_Read_Convec_EC( tmm, archivekey, tday, t1, t2, lli, levi, & gph, gph_tmi, & sp, & entu, entu_tmi, entd, entd_tmi, & detu, detu_tmi, detd, detd_tmi, & status ) use Phys , only : convec_mfldet_to_entdet use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(in) :: gph(:,:,:) ! m (half levels) type(TMeteoInfo), intent(in) :: gph_tmi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: entu(:,:,:), entd(:,:,:) ! kg/m2/s (full levels) type(TMeteoInfo), intent(out) :: entu_tmi, entd_tmi real, intent(out) :: detu(:,:,:), detd(:,:,:) ! kg/m2/s (full levels) type(TMeteoInfo), intent(out) :: detu_tmi, detd_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Convec_EC' ! --- local ------------------------------- integer :: lout real, allocatable :: ll(:,:,:) real, allocatable :: udmf(:,:,:) real, allocatable :: ddmf(:,:,:) real, allocatable :: uddr(:) real, allocatable :: dddr(:) integer :: i, j, k ! --- begin ------------------------------- call goLabel(rname) ! number of levels in output arrays: lout = size(entu,3) ! check ... if ( (size(entd,3)/=lout) .or. & (size(detu,3)/=lout) .or. (size(detd,3)/=lout) ) then write (gol,'("output arrays should have same number of levels:")'); call goErr write (gol,'(" entu : ",i4)') size(entu,3); call goErr write (gol,'(" entd : ",i4)') size(entd,3); call goErr write (gol,'(" detu : ",i4)') size(detu,3); call goErr write (gol,'(" detd : ",i4)') size(detd,3); call goErr TRACEBACK; status=1; return end if ! full level output: allocate( ll(lli%nlon,lli%nlat,levi%nlev ) ) ! generalized for lout=lmax_conv < levi%nlev allocate( udmf(lli%nlon,lli%nlat,lout+1) ) allocate( ddmf(lli%nlon,lli%nlat,lout+1) ) ! detrainment rates: allocate( uddr(lout) ) allocate( dddr(lout) ) #ifdef with_ec_aver ! ~~ archived fields are average over last hour ! updraught mass flux, half levels ! call ReadField( tmm, archivekey, 'MUMF', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entu_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... udmf(:,:,2:lout+1) = ll(:,:,1:lout) udmf(:,:,1) = 0. ! updraught detraintment rate call ReadField( tmm, archivekey, 'MUDR', 'kg/m3/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detu_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... detu = ll(:,:,1:lout) ! downdraught mass flux, half levels ! call ReadField( tmm, archivekey, 'MDMF', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entd_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... ddmf(:,:,2:lout+1) = ll(:,:,1:lout) ddmf(:,:,1) = 0. ! downdraught detraintment rate call ReadField( tmm, archivekey, 'MDDR', 'kg/m3/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detd_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... detd = ll(:,:,1:lout) #else ! ~~ fields are accumulated since start of run, ! deaccumulated on reading ! updraught mass flux, half levels ! call ReadField( tmm, archivekey, 'UDMF', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entu_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... udmf(:,:,2:lout+1) = ll(:,:,1:lout) udmf(:,:,1) = 0. ! updraught detraintment rate call ReadField( tmm, archivekey, 'UDDR', 'kg/m3/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detu_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... detu = ll(:,:,1:lout) ! downdraught mass flux, half levels ! call ReadField( tmm, archivekey, 'DDMF', 'kg/m2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, entd_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... ddmf(:,:,2:lout+1) = ll(:,:,1:lout) ddmf(:,:,1) = 0. ! downdraught detraintment rate call ReadField( tmm, archivekey, 'DDDR', 'kg/m3/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, detd_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... detd = ll(:,:,1:lout) #endif ! convert from flux/detr to entr/detr do j = 1, lli%nlat do i = 1, lli%nlon ! copy detrainment rates: uddr = detu(i,j,:) dddr = detd(i,j,:) ! convert: call convec_mfldet_to_entdet( 'u', lout, gph(i,j,1:lout+1), & udmf(i,j,1:lout+1), uddr , ddmf(i,j,1:lout+1), dddr , & entu(i,j,: ), detu(i,j,:), entd(i,j,: ), detd(i,j,:), status ) IF_NOTOK_RETURN(status=1) end do end do ! info ... call AddHistory( entu_tmi, 'converted mflux/detr to entr/detr', status ) call AddHistory( detu_tmi, 'converted mflux/detr to entr/detr', status ) call AddHistory( entd_tmi, 'converted mflux/detr to entr/detr', status ) call AddHistory( detd_tmi, 'converted mflux/detr to entr/detr', status ) ! clear deallocate( ll ) deallocate( udmf ) deallocate( ddmf ) deallocate( uddr ) deallocate( dddr ) ! ok call goLabel() status = 0 end subroutine tmm_Read_Convec_EC #endif ! ***************************************************************** ! ! Original parameters archived in MARS: ! ! References ! ERA-Interim archive plan ! http://www.ecmwf.int/publications/library/ecpublications/_pdf/era/era_report_series/RS_1_v2.pdf ! ! Parameter Surfaces Code Units ! ! turbulent diffusion coefficient for heat half levels 162.109 m2 (accum) ! subroutine tmm_Read_Diffus_EC( tmm, archivekey, tday, t1, t2, lli, levi, & sp, kzz, kzz_tmi, status ) use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli type(TLevelInfo), intent(in) :: levi real, intent(out) :: sp(:,:) ! Pa real, intent(out) :: kzz(:,:,:) ! m2/s (half levels) type(TMeteoInfo), intent(out) :: kzz_tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_Diffus_EC' ! --- local ------------------------------- integer :: lout real, allocatable :: ll(:,:,:) ! --- begin ------------------------------- ! number of full levels in output arrays: lout = size(kzz,3) - 1 ! full level output: allocate( ll(lli%nlon,lli%nlat,levi%nlev ) ) #ifdef with_ec_aver ! updraught mass flux, half levels, average over last hour: call ReadField( tmm, archivekey, 'MTDCH', 'm2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, kzz_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... kzz(:,:,2:lout+1) = ll(:,:,1:lout) kzz(:,:,1) = 0.0 #else ! updraught mass flux, half levels, accumulated from start of run: call ReadField( tmm, archivekey, 'K', 'm2/s', tday, t1, t2, & lli, 'n', levi, 'n', sp, ll, kzz_tmi, status ) IF_NOTOK_RETURN(status=1) ! store ... kzz(:,:,2:lout+1) = ll(:,:,1:lout) kzz(:,:,1) = 0.0 #endif ! info ... call AddHistory( kzz_tmi, 'implicit zero diffusion coefficient at surface', status ) ! clear deallocate( ll ) ! ok status = 0 end subroutine tmm_Read_Diffus_EC ! ================================================================== ! === ! === Olsson surface roughness ! === ! ================================================================== subroutine tmm_Read_SR_OLS( tmm, archivekey, tday, t1, t2, & lli, ll, tmi, status ) use GO , only : Get, goSplitLine use Grid , only : Init, Done, Interpol use tmm_info, only : TMeteoInfo, Init, AddHistory ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(:,:) type(TMeteoInfo), intent(out) :: tmi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_Read_SR_OLS' ! grid size integer, parameter :: nlon = 361, nlat = 181 ! --- local ------------------------------- character(len=10) :: sourcetype character(len=256) :: sourcename integer :: month character(len=256) :: fname logical :: exist, opened integer :: fu type(TllGridInfo) :: lli_ols real :: sr_ols(nlon,nlat) integer :: i, j ! --- begin ------------------------------- call goLabel(rname) ! split source key in type and name: call goSplitLine( archivekey, sourcetype, ':', sourcename, status ) IF_NOTOK_RETURN(status=1) ! input TMPP fields or raw prism fields ? select case ( sourcetype ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! standard ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'standard' ) ! dummy values: ll = 0.0 ! set history info call Init( tmi, 'srols', 'm', status ) call AddHistory( tmi, 'standard', status ) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! read directly from hdf file ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'tmpp', 'tm5-hdf', 'tm5-nc' ) ! read from surf file: call ReadField( tmm, archivekey, 'srols', 'm', tday, t1, t2, & lli, 'n', ll, tmi, status ) IF_NOTOK_RETURN(status=1) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! convert raw data: ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case ( 'ecmwf-tmpp', 'ecmwf-tm5', 'ncep-cdc', 'ncep-gfs' ) !call wrtgol( ' tmm read : srols ', t1, ' - ', t2 ); call goPr ! info ... call Init( tmi, 'srols', 'm', status ) ! month call Get( t1, month=month ) ! write filename: write (fname,'(a,"/SR_OLSSON_360_180_",i2.2,".d")') trim(tmm%input_dir), month ! info ... call AddHistory( tmi, 'file=='//trim(fname), status ) ! exist ? inquire( file=fname, exist=exist ) if ( .not. exist ) then write (fname,'(a,"/OLSSON-SR_OLSSON_360_180_",i2.2,".d")') trim(tmm%input_dir), month write(gol,*)" PLSPLS - using old filename for SR-OLSSON"; call goPr inquire( file=fname, exist=exist ) if ( .not. exist ) then write (gol,'("Olsson SR file not found:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return endif end if ! select free file unit: fu = 1234 do inquire( unit=fu, opened=opened ) if ( .not. opened ) exit fu = fu + 1 end do ! open data file: open( fu, file=trim(fname), form='formatted', status='old', iostat=status ) if (status/=0) then write (gol,'("while opening olsson data file:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! read field: read (fu,*,iostat=status) sr_ols if (status/=0) then write (gol,'("while reading from olsson data file:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! close file: close( fu, iostat=status ) if (status/=0) then write (gol,'("while closing olsson data file:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! setup grid definition: ! lon : -180.0 -179.0 .. 180.0 ( 1 deg resolution; 360 points; date line twice) ! lat : -90.0 -89.0 .. 90.0 ( 1 deg resolution; 180 points; includes poles) call Init( lli_ols, -180.00, 360.0/(nlon-1), nlon, -90.0, 180.0/(nlat-1), nlat, status ) IF_NOTOK_RETURN(status=1) ! interpol do j = 1, lli%nlat do i = 1, lli%nlon call Interpol( lli_ols, sr_ols, lli%lon_deg(i), lli%lat_deg(j), ll(i,j) ) end do end do ! info ... call AddHistory( tmi, 'horizontal_interpolation', status ) ! done call Done( lli_ols, status ) IF_NOTOK_RETURN(status=1) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! error ... ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ case default write (gol,'("unsupported source type `",a,"`")') trim(sourcetype); call goErr TRACEBACK; status=1; return end select ! ok call goLabel() status = 0 end subroutine tmm_Read_SR_OLS ! ! ================================================================== ! ! === ! ! === pv/theta -> eqv.lat. ! ! === ! ! ================================================================== ! ! ! subroutine tmm_ReadEqvLat( tmm, archivekey, & ! tday, t1, t2, & ! lli, levi, & ! sp, pv, theta, eqvlatb, eqvinds, & ! status ) ! ! use GO, only : TDate, operator(-), operator(+), operator(/), goSplitLine ! use Grid, only : TllGridInfo, TLevelInfo ! ! use tmm_mf , only : ReadEqvLatStuff ! use tmm_info, only : TMeteoInfo ! ! ! --- in/out -------------------------------- ! ! type(TTmMeteo), intent(inout) :: tmm ! ! character(len=*), intent(in) :: archivekey ! ! type(TDate), intent(in) :: tday, t1, t2 ! ! type(TllGridInfo), intent(in) :: lli ! type(TLevelInfo), intent(in) :: levi ! ! real, intent(out) :: sp(:,:) ! Pa ! real, intent(out) :: pv(:,:,:) ! real, intent(out) :: theta(:,:,:) ! real, intent(out) :: eqvlatb(:,:) ! integer, intent(out) :: eqvinds(:,:,:) ! ! integer, intent(out) :: status ! ! ! --- const -------------------------------------- ! ! character(len=*), parameter :: rname = mname//'/tmm_ReadEqvLat' ! ! ! --- local ------------------------------- ! ! character(len=10) :: sourcetype ! character(len=256) :: sourcename ! ! integer :: imf ! ! type(TMeteoInfo) :: tmi ! ! ! --- begin ------------------------------- ! ! ! split source key in type and name: ! call goSplitLine( archivekey, sourcetype, ':', sourcename ) ! ! ! input TMPP fields or raw prism fields ? ! select case ( sourcetype ) ! ! case ( 'tmpp' ) ! ! ! pv valid for [t1,t2] ! call ReadField( tmm, archivekey, 'PVo', tday, t1, t2, & ! lli, 'n', levi, 'n', sp, pv, tmi, status ) ! IF_NOTOK_RETURN(status=1) ! ! ! theta valid for [t1,t2] ! call ReadField( tmm, archivekey, 'theta', tday, t1, t2, & ! lli, 'n', levi, 'n', sp, theta, tmi, status ) ! IF_NOTOK_RETURN(status=1) ! ! ! ! ! eqv lat bounds and indices ! ! ! ! select (eventually retrieve first) the meteo file with this param: ! call SelectMF( tmm, 'i', archivekey, 'eqvlatb', tday, t1, t2, imf, status ) ! IF_NOTOK_RETURN(status=1) ! ! ! ! read from selected file ! call ReadEqvLatStuff( tmm%mf(imf), t1, t2, eqvlatb, eqvinds, status ) ! IF_NOTOK_RETURN(status=1) ! ! case default ! ! write (gol,'("unsupported source type `",a,"`")') trim(sourcetype) ! TRACEBACK; status=1; return ! ! end select ! ! ! ok ! status = 0 ! ! end subroutine tmm_ReadEqvLat ! ########################################################################################## ! ### ! ### output ! ### ! ########################################################################################## ! ! call WriteField( tmmd, 'od-fc-ml60-glb3x2', 'T', 'K', tday, t1, t2, & ! lli, 'n', sp, status ) ! subroutine tmm_WriteField_2d( tmm, archivekey, & tmi, paramkey, unit, tday, t1, t2, & lli, nuv, ll, status ) use tmm_mf , only : WriteRecord use tmm_info, only : TMeteoInfo ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TMeteoInfo), intent(in) :: tmi character(len=*), intent(in) :: paramkey character(len=*), intent(in) :: unit type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli character(len=1), intent(in) :: nuv real, intent(inout) :: ll(:,:) integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_WriteField_2d' ! --- local ---------------------------------------- integer :: imf ! --- begin ---------------------------------- !call wrtgol( 'tmm write : '//trim(paramkey)//' ', t1, ' - ', t2 ); call goPr ! check shape of grid: if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) ) then write (gol,'("2d array does not mach with grid definition:")'); call goErr write (gol,'(" param : ",a )') paramkey ; call goErr write (gol,'(" lli : ",i3," x ",i3)') lli%nlon, lli%nlat; call goErr write (gol,'(" nuv : ",a )') nuv; call goErr write (gol,'(" ll : ",i3," x ",i3)') shape(ll); call goErr TRACEBACK; status=1; return end if ! select index of already open meteo file or setup access to new one; call SelectMF( tmm, 'o', archivekey, paramkey, tday, t1, t2, imf, status ) IF_NOTOK_RETURN(status=1) ! write call WriteRecord( tmm%mf(imf), tmi, paramkey, unit, tday, t1, t2, & lli, nuv, ll, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine tmm_WriteField_2d ! ! call WriteField( tmmd, 'od-fc-ml60-glb3x2', 'T', 'K', tday, t1, t2, & ! lli, 'n', levi, spm, temper, status ) ! subroutine tmm_WriteField_3d( tmm, archivekey, & tmi, spname, paramkey, unit, tday, t1, t2, & lli, nuv, levi, nw, sp, ll, status )!, & !nlev ) use tmm_mf , only : WriteRecord use tmm_info, only : TMeteoInfo ! --- in/out -------------------------------- type(TTmMeteo), intent(inout) :: tmm character(len=*), intent(in) :: archivekey type(TMeteoInfo), intent(in) :: tmi character(len=*), intent(in) :: spname character(len=*), intent(in) :: paramkey character(len=*), intent(in) :: unit type(TDate), intent(in) :: tday, t1, t2 type(TllGridInfo), intent(in) :: lli character(len=1), intent(in) :: nuv type(TLevelInfo), intent(in) :: levi character(len=1), intent(in) :: nw real, intent(in) :: sp(:,:) ! Pa real, intent(inout) :: ll(:,:,:) integer, intent(out) :: status !integer, intent(in), optional :: nlev ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/tmm_WriteField_3d' ! --- local ---------------------------------------- integer :: imf ! --- begin ---------------------------------- !call wrtgol( 'tmm write : '//trim(paramkey)//' ', t1, ' - ', t2 ); call goPr ! check shape of grid: if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) .or. & ((nuv == 'n') .and. ((size(sp,1) /= lli%nlon ) .or. (size(sp,2) /= lli%nlat ))) .or. & ((nuv == 'u') .and. ((size(sp,1) /= lli%nlon+1) .or. (size(sp,2) /= lli%nlat ))) .or. & ((nuv == 'v') .and. ((size(sp,1) /= lli%nlon ) .or. (size(sp,2) /= lli%nlat+1))) .or. & ((nw == 'n') .and. (size(ll,3) > levi%nlev )) .or. & ((nw == 'w') .and. (size(ll,3) > levi%nlev+1)) ) then write (gol,'("3d arrays do not match with grid definition:")'); call goErr write (gol,'(" param : ",a )') paramkey; call goErr write (gol,'(" lli : ",i3," x ",i3 )') lli%nlon, lli%nlat; call goErr write (gol,'(" nuv : ",a )') nuv; call goErr write (gol,'(" levi : ",i3 )') levi%nlev; call goErr write (gol,'(" nw : ",a )') nw; call goErr write (gol,'(" sp : ",i3," x ",i3 )') shape(sp); call goErr write (gol,'(" ll : ",i3," x ",i3," x ",i3)') shape(ll); call goErr TRACEBACK; status=1; return end if ! select index of already open meteo file or setup access to new one; call SelectMF( tmm, 'o', archivekey, paramkey, tday, t1, t2, imf, status ) IF_NOTOK_RETURN(status=1) ! write call WriteRecord( tmm%mf(imf), tmi, spname, paramkey, unit, tday, t1, t2, & lli, nuv, levi, nw, sp, ll, status )!, & !nlev=nlev ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine tmm_WriteField_3d end module TMM