mo_diag_tendency.f90
MODULE mo_diag_tendency
! Authors:
5: ! I. Kirchner, MPI, October 1998 -
! add some fields for separation of terms in diagnostic equations
! ***************** Interface to ECHAM *************************
!
10: ! *ECHAM* *DIAGNOSE*
!
! CONTROL ----> DIAG_Init(1) general initialization
! +---> DIAG_Rerun('R') read rerun data
! STEPON -----> DIAG_Init(30) initialize transform buffer
15: ! +--> DIAG_Write store data in model output
! +---> DIAG_Check
! +---> DiagWriteSP
! +--> DIAG_Init(-1) clear diagnostic memory
! +---> DIAG_Rerun('W') store rerun data
20: ! SCAN1SL ----> DIAG_fftd
! +---> DIAG_Init(40) initialise buffer
! +--> DIAG_sym1
! +--> DIAG_ltd
!
25: ! count terms in additional arrays
! *spectral* *Gaussian*
!
! DYN adiabatic terms
! TF2 time filter
30: ! |
! TF1 time filter
! GPC
! +--> PHYSC diabatic terms
! | +---> M_RADHEAT long/short wave
35: ! +--> SI1 semi-implicit terms
! DIAG_fftd
! SI2 semi-implicit terms
! DIAG_sym1
! DIAG_ltd
40: ! SCCD semi-implicit terms
! SCCTP semi-implicit terms
! HDIFF horizontal diffusion
!
!**************************************************************
45:
USE mo_doctor, ONLY: nout
USE mo_exception, ONLY: finish
USE mo_start_dataset, ONLY: lres, nstep
USE mo_year, ONLY: ic2ymd
50:
USE mo_truncation, ONLY: nmp, nnp
USE mo_control, ONLY: &
& nlev, nlevp1, nsp, n2sp, nlp2, ngl, nmp1, nhgl, &
& lptime, twodt, nrow, nlon, nptime, nvclev, vct
55:
USE mo_fft
USE mo_legendre, ONLY: pnmd, anmd, rnmd
USE mo_post, ONLY: nunitsp, nunitdf
60: USE mo_grib
USE mo_memory_sp, ONLY: sd, svo, stp
IMPLICIT NONE
65: CHARACTER :: diarev*20 = ' Rev. 2.3 10-OCT-99' ! label for modul revision
! LTDIAG additional switch in &RUNCTL
! .true. ==> addional diagnostics
!
70: ! DIAG_Init allocate/deallocate memory, initialize diagnostic arrays
! DIAG_fftd Fourier transformation of diagnostics arrays
! DIAG_sym1 separation in symetric/antisymetric part
! DIAG_ltd legendre transform of diagnostic arrays
! DIAG_Write store diagnostic fields
75: ! DIAG_Rerun rerun data handling
! special diagnostic fields
REAL, ALLOCATABLE :: & ! ARRAYS in spectral space
80: & pdvor (:,:,:,:), & ! terms for vorticity equation
& pddiv (:,:,:,:), & ! terms for divergence equation
& pdtem (:,:,:,:), & ! terms for temperature equation
& pdprs (:,:,:), & ! terms for surface pressure equation
& pdprl (:,:,:) ! surface pressure term vertical integration
85:
REAL, ALLOCATABLE :: & ! old value for accumulation check
& pdovor (:,:,:), pdodiv (:,:,:), pdotep (:,:,:)
LOGICAL :: ltdiag2 ! store state of diagnostics
90: INTEGER, PARAMETER :: & ! number of terms in equation
& ndvor=8, nddiv=8, ndtem=13, ndprs=3
LOGICAL :: ldinit ! state of initialization
REAL, ALLOCATABLE :: & ! ARRAYS in grid point space
95: & ptfh1 (:,:,:,:), & ! time integration parts
& ptfh2 (:,:)
REAL, ALLOCATABLE :: & ! ARRAYS for counting
& pdiga (:,:,:,:), zdiga (:), pdigaa(:,:,:,:,:), pdigas(:,:,:,:,:), &
100: & pdsga (:,:,:), zdsga (:), pdsgaa(:,:,:,:), pdsgas(:,:,:,:), &
& pdigb (:,:,:,:), pdigba(:,:,:,:,:), pdigbs(:,:,:,:,:), &
& pdigs (:,:,:,:), zdigs (:), pdigsa(:,:,:,:,:), pdigss(:,:,:,:,:), &
& pdsgs (:,:), zdsgs (:), pdsgsa(:,:,:), pdsgss(:,:,:)
105:
INTEGER,PARAMETER :: ldunit=90 ! unit for rerun tendency data
! insert the new codes here, if GRIB1 works fine into the postprecessing
!
110: ! codes for diagnostic terms
! vorticity
INTEGER, PARAMETER :: ncvor(ndvor)=(/41,42,43,44,45,46,47,48/)
!integer, parameter :: ncvor(ndvor)=(/155,156,157,158,159,160,161/)
! 1 041/155 dynamic tendencies horizontal part
115: ! 2 042/156 vertical advection
! 3 043/157 vertical dissipation due to VDIFF
! 4 044/158 vertical dissipation due to GWDRAG
! 5 045/159 vertical dissipation due to CUCALL
! 6 046/160 timefilter
120: ! 7 047/161 semi-implicit time integration part
! 8 048/162 horizontal diffusion
!
! divergence
INTEGER, PARAMETER :: ncdiv(nddiv)=(/61,62,63,64,65,66,67,68/)
125: !integer, parameter :: ncdiv(nddiv)=(/170,171,172,173,174,175,176,177/)
! 1 061/170 dynamic tendencies horizontal part
! 2 062/171 vertical advection
! 3 063/172 vertical dissipation due to VDIFF
! 4 064/173 vertical dissipation due to GWDRAG
130: ! 5 065/174 vertical dissipation due to CUCALL
! 6 066/175 timefilter
! 7 067/176 semi-implicit time integration part
! 8 068/177 horizontal diffusion
!
135: ! temperature
INTEGER, PARAMETER :: nctem(ndtem)=(/81,82,83,84,85,86,87,88,89,90,91,92,93/)
!integer, parameter :: nctem(ndtem)=(/185,186,187,188,189,190,191,192,193,194,195,196,197/)
! 1 081/185 dynamic tendencies horizontal part
! 2 082/186 vertical advection
140: ! 3 083/187 energy conversion
! 4 084/188 radiation processes
! 5 085/189 vertical dissipation due to VDIFF
! 6 086/190 vertical dissipation due to GWDRAG
! 7 087/191 cond, prec, eva due to CUCALL
145: ! 8 088/192 cond, prec, eva due to CUCOND
! 9 089/193 timefilter
! 10 090/194 semi-implicit time integration part
! 11 091/195 horizontal diffusion
! 12 092/196 longwave radiation
150: ! 13 093/197 shortwave radiation
!
! surface pressure
INTEGER, PARAMETER :: ncprs(ndprs)=(/101,102,103/)
!integer, parameter :: ncprs(ndprs)=(/201,202,203/)
155: ! 1 101/201 divergence term
! 2 102/202 timefilter
! 3 103/203 semi-implicit time integration part
INTEGER, PARAMETER :: ncprl = 100
!integer, parameter :: ncprl = 200
160: ! 1 100/200 horizontal divergence in the layers
CONTAINS
!======================================================================
165: SUBROUTINE DIAG_Init(itype)
INTEGER, INTENT(in) :: itype ! function separator
SELECT CASE(itype)
170: CASE(1) ! startup initialization
ldinit = .TRUE.
ltdiag2 = .FALSE.
175: WRITE (nout,*) &
& ' Diagnostics of tendency terms ',TRIM(diarev),' (kirchner@dkrz.de)'
! -----------------------------------------------------------
! allocate spectral space arrays
!
180: ! diagnostic arrays all accumulated
!
! g .... accumulated in grid point space
! gs ... parts from gridpoint space, accumulated in spectral space
! s .... accumulated in spectral space
185: !
! spectral terms
!
! PDVOR vorticity equation
ALLOCATE (pdvor (nlev,2,nsp,ndvor)); pdvor(:,:,:,:) = 0.0
190: ALLOCATE (pdovor (nlev,2,nsp))
!
! g 1 horizontal advection, pressure gradient, coriolisterm (DYN)
! g 2 vertical advection (DYN)
! g 3 vertical diffusion due to impuls (VDIFF)
195: ! g 4 gravity wave drag (GWDRAG)
! g 5 moisture mass flux (CUCALL)
! g 6 timefilter
! gs 7 semi-implicit part of time integration
! s 8 horizontal diffusion
200: !
! PDDIV divergence equation
ALLOCATE (pddiv (nlev,2,nsp,nddiv)); pddiv(:,:,:,:) = 0.0
ALLOCATE (pdodiv (nlev,2,nsp))
!
205: ! g 1 horizontal advection, coriolisterm, pressure gradient, G-term
! g 2 vertical advection
! g 3 vertical diffusion due to imuls (VDIFF)
! g 4 gravity wave drag (GWDRAG)
! g 5 moisture mass flux (CUCALL)
210: ! g 6 timefilter
! gs 7 semi-implicit part of time integration
! s 8 horizontal diffusion
!
! PDTEM temperature
215: ALLOCATE (pdtem (nlev,2,nsp,ndtem)); pdtem(:,:,:,:) = 0.0
ALLOCATE (pdotep (nlevp1,2,nsp))
!
! g 1 horizontal advection
! g 2 vertical advection
220: ! g 3 energy conversion
! g 4 radiation (RADHEAT)
! g 5 vertical diffusion due to turbulence (VDIFF)
! g 6 gravity wave drag (GWDRAG)
! g 7 convection (CUCALL)
225: ! g 8 large scale cloud processes (COND)
! g 9 timefilter
! gs 10 semi-implicit part of time integration
! s 11 horizontal diffusion
! g 12 longwave radiation
230: ! g 13 shortwave radiation
!
! PDPRL divergence for each layer
ALLOCATE (pdprl (nlev,2,nsp)); pdprl(:,:,:) = 0.0
! g convergence in each layer
235: !
! PDPRS pressure equation
ALLOCATE (pdprs (2,nsp,ndprs)); pdprs(:,:,:) = 0.0
!
! g 1 vertical integrated convergence
240: ! g 2 timefilter
! s 3 semi-implicit part of time integration
!
! -----------------------------------------------------------
! allocate gridpoint space arrays
245: !
! PTFH1 memory of timefilter
! updated in TF1
ALLOCATE (ptfh1 (nlp2,nlev,3,ngl)); ptfh1(:,:,:,:) = 0.0
! 1 ... vorticity
250: ! 2 ... divergence
! 3 ... temperature
! PTFH2 pressure
! updated in TF1
ALLOCATE (ptfh2 (nlp2, ngl)); ptfh2(:,:) = 0.0
255: !
! -----------------------------------------------------------
! workspace for accumulation of terms in grid point space
!
! Index
260: ! ............. vorticity and divergenc equation
! VOM = Fu, VOL = Fv
!
! Fu (VOM parts in GPC)
! DYN 1 coriolisterm and pressure tendency of Fu
265: ! DYN 2 vertical advection of Fu
! PHYSC 3 diffusion due to impuls Pu VDIFF
! PHYSC 4 diffusion due to gravity wave drag Pu GWDRAG
! PHYSC 5 diffusion due to mass flux Pu CUCALL
!
270: ! Fv (VOL parts in GPC)
! DYN 6 coriolisterm and pressure tendency of Fv
! DYN 7 vertical advection of Fv
! PHYSC 8 diffusion due to impuls Pv VDIFF
! PHYSC 9 diffusion due to gravity wave drag Pv GWDRAG
275: ! PHYSC 10 diffusion due to mass flux Pv CUCALL
!
! DYN 11 potential and kinetic energy term G
!
! ............. temperature tendency equation
280: ! DYN 12 horizontal advection term
! DYN 13 vertical advection term
! DYN 14 energy conversion
! PHYSC 15 RADHEAT radiation tendency
! PHYSC 16 VDIFF turbulence
285: ! PHYSC 17 GWDRAG gravity wave drag
! PHYSC 18 CUCALL convective condensation
! PHYSC 19 COND large+scale condensation
!
! ............. pressure tendency equation
290: ! DYN 20 level dependend divergence part,
!
! TF2 21 timefilter of vorticity
! TF2 22 timefilter of divergence
! TF2 23 timefilter of temperature
295: !
! RADHEAT 24 longwave radiation
! RADHEAT 25 shortwave radiation
!
ALLOCATE (pdiga (nlp2,nlev,25,ngl)); pdiga(:,:,:,:) = 0.0
300: !
ALLOCATE (zdiga (nlp2*nlev*25)) ! workspace fourier transform
ALLOCATE (pdigaa (nlev,2,nmp1,25,nhgl))
ALLOCATE (pdigas (nlev,2,nmp1,25,nhgl))
!
305: ! one level arrays of surface pressure
! 1 last level total integral
! 2 timefilter for pressure
!
ALLOCATE (pdsga (nlp2,2,ngl)); pdsga(:,:,:) = 0.0
310: !
ALLOCATE (zdsga (nlp2*2)) ! workspace fourier transform
ALLOCATE (pdsgaa (2,nmp1,2,nhgl))
ALLOCATE (pdsgas (2,nmp1,2,nhgl))
!
315: ! array for d/dlambda derivatives calculated in fourierspace
! corresponding to pdiga(*,1...10,*)
!
ALLOCATE (pdigb (nlp2,nlev,10,ngl))
ALLOCATE (pdigba (nlev,2,nmp1,10,nhgl))
320: ALLOCATE (pdigbs (nlev,2,nmp1,10,nhgl))
!
! local memory buffer for accumulation of semi-implicit parts
! fraction solved in grid point space
!
325: ! 1 vorticity and implicit and explicit part (L)
! 2 divergence
! 3 temperatur and pressure
! 4 implicit part of vorticity (M) used in Fourierspace
ALLOCATE (pdigs (nlp2,nlev,4,ngl))
330: !
ALLOCATE (zdigs (nlp2*nlev*3)) ! workspace fourier transform
ALLOCATE (pdigsa (nlev,2,nmp1,4,nhgl))
ALLOCATE (pdigss (nlev,2,nmp1,4,nhgl))
!
335: ALLOCATE (pdsgs (nlp2,ngl))
!
ALLOCATE (zdsgs (nlp2)) ! workspace fourier transform
ALLOCATE (pdsgsa (2,nmp1,nhgl))
ALLOCATE (pdsgss (2,nmp1,nhgl))
340: !
! special initialization for restart
IF (lres) CALL DIAG_Rerun('R') ! read old rerun data
CASE(-1) ! get memory free
345:
CALL DIAG_Rerun('W') ! store rerun data
DEALLOCATE (pdvor)
DEALLOCATE (pddiv)
350: DEALLOCATE (pdtem)
DEALLOCATE (pdprs)
DEALLOCATE (pdprl)
DEALLOCATE (pdovor)
355: DEALLOCATE (pdodiv)
DEALLOCATE (pdotep)
DEALLOCATE (ptfh1)
DEALLOCATE (ptfh2)
360:
DEALLOCATE (pdiga)
DEALLOCATE (zdiga)
DEALLOCATE (pdigaa)
DEALLOCATE (pdigas)
365:
DEALLOCATE (pdsga)
DEALLOCATE (zdsga)
DEALLOCATE (pdsgaa)
DEALLOCATE (pdsgas)
370:
DEALLOCATE (pdigb)
DEALLOCATE (pdigba)
DEALLOCATE (pdigbs)
375: DEALLOCATE (pdigs)
DEALLOCATE (zdigs)
DEALLOCATE (pdigsa)
DEALLOCATE (pdigss)
380: DEALLOCATE (pdsgs)
DEALLOCATE (zdsgs)
DEALLOCATE (pdsgsa)
DEALLOCATE (pdsgss)
385: CASE(30) ! initialize legendre buffer in SCAN1SL
! necessary at every time step
pdigaa(:,:,:,:,:) = 0.0
pdigas(:,:,:,:,:) = 0.0
390: pdsgaa(:,:,:,:) = 0.0
pdsgas(:,:,:,:) = 0.0
pdigb (:,:,:,:) = 0.0
pdigba(:,:,:,:,:) = 0.0
395: pdigbs(:,:,:,:,:) = 0.0
pdigs (:,:,:,:) = 0.0
pdigsa(:,:,:,:,:) = 0.0
pdigss(:,:,:,:,:) = 0.0
400:
pdsgs (:,:) = 0.0
pdsgsa(:,:,:) = 0.0
pdsgss(:,:,:) = 0.0
405: CASE(40) ! initialise local fourierbuffer in DIAG_fft
! necessary at every latitude
zdiga(:) = 0.0
zdsga(:) = 0.0
410: zdigs(:) = 0.0
zdsgs(:) = 0.0
CASE default
WRITE (nout,*) ' DIAG_Init: type not implemented, ITYPE= ',itype
415: END SELECT
RETURN
END SUBROUTINE DIAG_Init
420: !======================================================================
SUBROUTINE DIAG_fftd ! -->>> insert in SCAN1SL after CALL FFTD
! transform the diagnostic arrays into fourierspace
INTEGER :: irow ! actual latitude index
425: INTEGER :: inc, isign
REAL :: zwork(nlp2*nlev*25) ! work space for transformation
#if (defined CRAY) || (defined SX)
INTEGER :: iamount
EXTERNAL util_reshape
430: #endif
INTRINSIC RESHAPE
irow = nrow(1)
435: CALL DIAG_init(40)
! Load input arrays
#if (defined CRAY) || (defined SX)
iamount = nlp2*nlev*3
440: CALL util_reshape(zdigs,pdigs(1,1,1,irow),iamount)
iamount = nlp2
CALL util_reshape(zdsgs,pdsgs(1 ,irow),iamount)
#else
zdigs(:) = RESHAPE(pdigs(:,:,1:3 ,irow),(/nlp2*nlev*3/))
445: zdsgs(:) = RESHAPE(pdsgs(: ,irow),(/nlp2/))
#endif
inc = 1
isign = -1
450:
CALL fft991cy(zdigs,zwork,trig,nfax,inc,nlp2,nlon,3*nlev,isign)
CALL fft991cy(zdsgs,zwork,trig,nfax,inc,nlp2,nlon,1,isign)
! reshape the transformed arrays
455: #if (defined CRAY) || (defined SX)
iamount = nlp2*nlev*3
CALL util_reshape(pdigs(1,1,1,irow),zdigs,iamount)
iamount = nlp2
CALL util_reshape(pdsgs(1 ,irow),zdsgs,iamount)
460: #else
pdigs(:,:,1:3 ,irow) = RESHAPE(zdigs(:),(/nlp2,nlev,3/))
pdsgs(: ,irow) = RESHAPE(zdsgs(:),(/nlp2/))
#endif
465: IF (lptime) THEN ! other terms during output time step
! Load input arrays
#if (defined CRAY) || (defined SX)
iamount = nlp2*nlev*25
470: CALL util_reshape(zdiga,pdiga(1,1,1,irow),iamount)
iamount = nlp2*2
CALL util_reshape(zdsga,pdsga(1,1 ,irow),iamount)
#else
zdiga(:) = RESHAPE(pdiga(:,:,:,irow),(/nlp2*nlev*25/))
475: zdsga(:) = RESHAPE(pdsga(:,: ,irow),(/nlp2*2/))
#endif
CALL fft991cy(zdiga,zwork,trig,nfax,inc,nlp2,nlon,25*nlev,isign)
CALL fft991cy(zdsga,zwork,trig,nfax,inc,nlp2,nlon,2,isign)
480:
! reshape the transformed arrays
#if (defined CRAY) || (defined SX)
iamount = nlp2*nlev*25
CALL util_reshape(pdiga(1,1,1,irow),zdiga,iamount)
485: iamount = nlp2*2
CALL util_reshape(pdsga(1,1 ,irow),zdsga,iamount)
#else
pdiga(:,:,:,irow) = RESHAPE(zdiga(:),(/nlp2,nlev,25/))
pdsga(:,: ,irow) = RESHAPE(zdsga(:),(/nlp2,2/))
490: #endif
ENDIF
RETURN
495: END SUBROUTINE DIAG_fftd
!======================================================================
SUBROUTINE DIAG_sym1 ! -->>> insert in SCAN1SL after CALL SYM1
! separate the diagnostics arrays
500: INTEGER :: ihrow ! half latitude index
INTEGER :: jl, jm, irow
irow = nrow(1)
ihrow = (irow+1)/2
505:
! even and odd components
IF (MOD(irow,2) /= 0) THEN ! northern hemisphere
DO jm = 1,nmp1
510: DO jl = 1,nlev
pdigss(jl,:,jm,:,ihrow) = pdigs(2*jm-1:2*jm,jl,:,irow)
pdigsa(jl,:,jm,:,ihrow) = pdigs(2*jm-1:2*jm,jl,:,irow)
ENDDO
pdsgss(:,jm,ihrow) = pdsgs(2*jm-1:2*jm,irow)
515: pdsgsa(:,jm,ihrow) = pdsgs(2*jm-1:2*jm,irow)
ENDDO
IF (lptime) THEN
DO jm = 1,nmp1
DO jl = 1,nlev
520: pdigas(jl,:,jm,:,ihrow) = pdiga(2*jm-1:2*jm,jl,:,irow)
pdigaa(jl,:,jm,:,ihrow) = pdiga(2*jm-1:2*jm,jl,:,irow)
pdigbs(jl,:,jm,:,ihrow) = pdigb(2*jm-1:2*jm,jl,:,irow)
pdigba(jl,:,jm,:,ihrow) = pdigb(2*jm-1:2*jm,jl,:,irow)
ENDDO
525: pdsgas(:,jm,:,ihrow) = pdsga(2*jm-1:2*jm,:,irow)
pdsgaa(:,jm,:,ihrow) = pdsga(2*jm-1:2*jm,:,irow)
ENDDO
ENDIF
ELSE ! southern hemisphere
530: DO jm = 1,nmp1
DO jl = 1,nlev
pdigss(jl,:,jm,:,ihrow) = 0.5*(pdigss(jl,:,jm,:,ihrow) + pdigs(2*jm-1:2*jm,jl,:,irow))
pdigsa(jl,:,jm,:,ihrow) = 0.5*(pdigsa(jl,:,jm,:,ihrow) - pdigs(2*jm-1:2*jm,jl,:,irow))
ENDDO
535: pdsgss(:,jm,ihrow) = 0.5*(pdsgss(:,jm,ihrow) + pdsgs(2*jm-1:2*jm,irow))
pdsgsa(:,jm,ihrow) = 0.5*(pdsgsa(:,jm,ihrow) - pdsgs(2*jm-1:2*jm,irow))
ENDDO
IF (lptime) THEN
DO jm = 1,nmp1
540: DO jl = 1,nlev
pdigas(jl,:,jm,:,ihrow) = 0.5*(pdigas(jl,:,jm,:,ihrow) + pdiga(2*jm-1:2*jm,jl,:,irow))
pdigaa(jl,:,jm,:,ihrow) = 0.5*(pdigaa(jl,:,jm,:,ihrow) - pdiga(2*jm-1:2*jm,jl,:,irow))
pdigbs(jl,:,jm,:,ihrow) = 0.5*(pdigbs(jl,:,jm,:,ihrow) + pdigb(2*jm-1:2*jm,jl,:,irow))
pdigba(jl,:,jm,:,ihrow) = 0.5*(pdigba(jl,:,jm,:,ihrow) - pdigb(2*jm-1:2*jm,jl,:,irow))
545: ENDDO
pdsgas(:,jm,:,ihrow) = 0.5*(pdsgas(:,jm,:,ihrow) + pdsga(2*jm-1:2*jm,:,irow))
pdsgaa(:,jm,:,ihrow) = 0.5*(pdsgaa(:,jm,:,ihrow) - pdsga(2*jm-1:2*jm,:,irow))
ENDDO
ENDIF
550: ENDIF
RETURN
END SUBROUTINE DIAG_sym1
!======================================================================
555: SUBROUTINE DIAG_ltd ! -->>> insert in SCAN1SL after CALL LTD
! perform legendre transform for diagnostic arrays
INTEGER :: ihrow ! half latitude index
INTEGER :: jm, ims, ins, is, jn, jh, iu
560:
ihrow = nrow(1)
DO jh = 1,2 ! 1: north, 2:south
iu = 2-jh
565: DO jm = 1,nmp1
ims = nmp(jm)-iu
ins = nnp(jm)+iu
DO jn=2,ins,2
is = ims+jn
570: IF (jh == 1) THEN ! calculations for northern hemisphere
!
! semi-implicit parts of vorticity
pdvor(:,:,is, 7) = pdvor(:,:,is, 7) + pdigss(:,:,jm,1,ihrow)*pnmd(is) &
& + pdigsa(:,:,jm,4,ihrow)*anmd(is)
575:
! explicit part of divergence, temperature and pressure
pddiv(:,:,is, 7) = pddiv(:,:,is, 7) + pdigss(:,:,jm,2,ihrow)*rnmd(is)
pdtem(:,:,is,10) = pdtem(:,:,is,10) + pdigss(:,:,jm,3,ihrow)*pnmd(is)
pdprs (:,is, 3) = pdprs (:,is, 3) + pdsgss (:,jm ,ihrow)*pnmd(is)
580: !
IF (lptime) THEN
! dynamic and physical tendencies
! vorticity
pdvor(:,:,is,1:5) = pdvor (:,:,is,1:5) + pdigbs(:,:,jm,6:10,ihrow)*pnmd(is) &
585: & + pdigaa(:,:,jm,1:5 ,ihrow)*anmd(is)
pdvor(:,:,is, 6) = pdvor (:,:,is, 6) + pdigas(:,:,jm,21 ,ihrow)*pnmd(is)
!
! divergence
pddiv(:,:,is,1:5) = pddiv (:,:,is,1:5) + pdigbs(:,:,jm,1:5 ,ihrow)*pnmd(is) &
590: & - pdigaa(:,:,jm,6:10,ihrow)*anmd(is)
pddiv(:,:,is, 6) = pddiv (:,:,is, 6) + pdigas(:,:,jm,22 ,ihrow)*pnmd(is)
! laplacian operation with G-term
pddiv(:,:,is, 1) = pddiv (:,:,is, 1) + pdigas(:,:,jm,11 ,ihrow)*rnmd(is)
!
595: ! temperature
pdtem(:,:,is, 1: 8) = pdtem(:,:,is, 1: 8) + pdigas(:,:,jm,12:19,ihrow)*pnmd(is)
pdtem(:,:,is, 9) = pdtem(:,:,is, 9) + pdigas(:,:,jm,23 ,ihrow)*pnmd(is)
pdtem(:,:,is,12:13) = pdtem(:,:,is,12:13) + pdigas(:,:,jm,24:25,ihrow)*pnmd(is)
!
600: ! pressure
pdprl(:,:,is) = pdprl(:,:,is) + pdigas(:,:,jm,20,ihrow)*pnmd(is)
pdprs (:,is,1:2) = pdprs (:,is,1:2) + pdsgas (:,jm, :,ihrow)*pnmd(is)
ENDIF
!
605: ELSE ! calculations for southern hemisphere
!
! semi-implicit parts of vorticity
pdvor(:,:,is, 7) = pdvor (:,:,is, 7) + pdigsa(:,:,jm,1,ihrow)*pnmd(is) &
& + pdigss(:,:,jm,4,ihrow)*anmd(is)
610: ! (-) hier ?? oder (+)
! explicit part divergence, temperature and pressure
pddiv(:,:,is, 7) = pddiv(:,:,is, 7) + pdigsa(:,:,jm,2,ihrow)*rnmd(is)
pdtem(:,:,is,10) = pdtem(:,:,is,10) + pdigsa(:,:,jm,3,ihrow)*pnmd(is)
615: pdprs (:,is, 3) = pdprs (:,is, 3) + pdsgsa (:,jm ,ihrow)*pnmd(is)
!
IF (lptime) THEN
! dynamic and physical tendencies
! vorticity
620: pdvor(:,:,is,1:5) = pdvor(:,:,is,1:5) + pdigba(:,:,jm,6:10,ihrow)*pnmd(is) &
& + pdigas(:,:,jm,1:5 ,ihrow)*anmd(is)
pdvor(:,:,is, 6) = pdvor(:,:,is, 6) + pdigaa(:,:,jm,21 ,ihrow)*pnmd(is)
!
! divergence
625: pddiv(:,:,is,1:5) = pddiv(:,:,is,1:5) + pdigba(:,:,jm,1:5 ,ihrow)*pnmd(is) &
& - pdigas(:,:,jm,6:10,ihrow)*anmd(is)
pddiv(:,:,is, 6) = pddiv(:,:,is, 6) + pdigaa(:,:,jm,22 ,ihrow)*pnmd(is)
pddiv(:,:,is, 1) = pddiv(:,:,is, 1) + pdigaa(:,:,jm,11 ,ihrow)*rnmd(is)
!
630: ! temperature
pdtem(:,:,is, 1: 8) = pdtem(:,:,is, 1: 8) + pdigaa(:,:,jm,12:19,ihrow)*pnmd(is)
pdtem(:,:,is, 9) = pdtem(:,:,is, 9) + pdigaa(:,:,jm,23 ,ihrow)*pnmd(is)
pdtem(:,:,is,12:13) = pdtem(:,:,is,12:13) + pdigaa(:,:,jm,24:25,ihrow)*pnmd(is)
!
635: ! pressure
pdprl(:,:,is) = pdprl(:,:,is) + pdigaa(:,:,jm,20,ihrow)*pnmd(is)
pdprs (:,is,1:2) = pdprs (:,is,1:2) + pdsgaa (:,jm, :,ihrow)*pnmd(is)
ENDIF
!
640: ENDIF
ENDDO
ENDDO
ENDDO
645: RETURN
END SUBROUTINE DIAG_ltd
!======================================================================
SUBROUTINE DIAG_Write
REAL :: tfact
650: INTEGER :: iterm, jlev, ierr, iret
REAL :: worksp(2,nsp)
EXTERNAL codegb5, pbwrite
655: ! proportionality factor
! all tendencies are stored as VALUE/sec
tfact = 1.0/(twodt*nptime)
! Finish definition GRIB block 1
660:
ksec1(1) = nudging_table
CALL set_output_time
665: ksec1(10) = year
ksec1(11) = month
ksec1(12) = day
ksec1(13) = hour
ksec1(14) = minute
670: ksec1(21) = century
ksec4(1) = n2sp
WRITE (nout,*) 'DIAG_Write header: century year mon day hour min ', &
675: & century,year,month,day,hour,minute
level_type = 109
ksec1(7) = level_type
680: ! perform a internal check
IF (ltdiag2) CALL DIAG_Check
DO iterm = 1,ndvor ! 1. vorticity equation
ksec1(6) = ncvor(iterm)
685: DO jlev = 1 , nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = pdvor(jlev,:,:,iterm)*tfact
CALL DiagWriteSP
690: ENDDO
ENDDO
pdvor(:,:,:,:) = 0.0
pdovor(:,:,:) = svo(:,:,:)
695: DO iterm = 1,nddiv ! 2. divergence equation
ksec1(6) = ncdiv(iterm)
DO jlev = 1 , nlev
level_p2 = jlev
ksec1(9) = level_p2
700: worksp(:,:) = pddiv(jlev,:,:,iterm)*tfact
CALL DiagWriteSP
ENDDO
ENDDO
pddiv(:,:,:,:) = 0.0
705: pdodiv(:,:,:) = sd(:,:,:)
DO iterm = 1,ndtem ! 3. temperature equation
ksec1(6) = nctem(iterm)
DO jlev = 1 , nlev
710: level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = pdtem(jlev,:,:,iterm)*tfact
CALL DiagWriteSP
ENDDO
715: ENDDO
pdtem(:,:,:,:) = 0.0
pdotep(:,:,:) = stp(:,:,:)
ltdiag2 = .TRUE. ! now diagnostic fields initialised
720:
ksec1(6) = ncprl ! 4. pressure equation
DO jlev = 1 , nlev
level_p2 = jlev
ksec1(9) = level_p2
725: worksp(:,:) = pdprl(jlev,:,:)*tfact
CALL DiagWriteSP
ENDDO
pdprl(:,:,:) = 0.0
730: level_type = 1 ! surface data
ksec1(7) = level_type
level_p2 = 0
ksec1(9) = level_p2
DO iterm = 1,ndprs
735: ksec1(6) = ncprs(iterm)
worksp(:,:) = pdprs(:,:,iterm)*tfact
CALL DiagWriteSP
ENDDO
pdprs(:,:,:) = 0.0
740:
! reset accumulated grid point arrays
pdiga(:,:,:,:) = 0.0
pdsga(:,:,:) = 0.0
745: RETURN
CONTAINS
SUBROUTINE DiagWriteSP
750: #ifdef EMOS
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
& ksec4, worksp, n2sp, kgrib, kleng, kword, 'C', ierr)
#else
CALL codegb5(worksp,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nvclev,&
755: & kgrib,klengo,kword,0,ierr)
#endif
CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
IF (iret /= kword*iobyte) &
& CALL finish('DIAG_Write','I/O error on output - disk full?')
760: RETURN
END SUBROUTINE DiagWriteSP
END SUBROUTINE DIAG_Write
!======================================================================
765: SUBROUTINE DIAG_Rerun(ctyp)
INTEGER :: ilevels, ilons2, ilats, insp, jl, flag
CHARACTER :: ctyp*1
770: REAL, ALLOCATABLE :: rvct(:)
SELECT CASE(ctyp)
CASE('r','R') ! read rerun data from file
REWIND(ldunit,err=100)
775:
READ(ldunit,err=100) ilevels, ilons2, ilats, insp
IF (ilevels /= nlev) THEN
WRITE (nout,*) 'DIAG_Rerun Warning: number of levels inconsistent'
780: ELSE IF (ilons2 /= nlp2) THEN
WRITE (nout,*) 'DIAG_Rerun Warning: number of longitudes inconsistent'
ELSE IF (ilats /= ngl) THEN
WRITE (nout,*) 'DIAG_Rerun Warning: number of latitudes inconsistent'
ELSE IF (insp /= nsp) THEN
785: WRITE (nout,*) 'DIAG_Rerun Warning: number of spectral coefficients inconsistent'
ELSE
! get coordiante system from file
ALLOCATE (rvct(2*nlev+2))
READ(ldunit,err=100) rvct
790: flag = 0
DO jl = 1,2*nlev+2
IF (rvct(jl) /= vct(jl)) flag = flag + 1
ENDDO
DEALLOCATE (rvct)
795: IF (flag /= 0) THEN
WRITE (nout,*) 'DIAG_Rerun Warning: inconsistent vertical model structure'
ELSE
READ(ldunit) pdvor
READ(ldunit) pddiv
800: READ(ldunit) pdtem
READ(ldunit) pdprl
READ(ldunit) pdprs
READ(ldunit) ptfh1
READ(ldunit) ptfh2
805: ldinit = .FALSE. ! no additional initialization for timefilter
WRITE (nout,*) 'DIAG_Rerun ... read restart data'
ENDIF
ENDIF
810: CASE('w','W') ! write restart data to file
REWIND(ldunit)
WRITE(ldunit) nlev,nlp2,ngl,nsp
WRITE(ldunit) vct(1:2*nlev+2)
WRITE(ldunit) pdvor
815: WRITE(ldunit) pddiv
WRITE(ldunit) pdtem
WRITE(ldunit) pdprl
WRITE(ldunit) pdprs
WRITE(ldunit) ptfh1
820: WRITE(ldunit) ptfh2
WRITE (nout,*) 'DIAG_Rerun ... wrote restart data'
END SELECT
825:
100 CONTINUE
RETURN
END SUBROUTINE DIAG_Rerun
830:
SUBROUTINE DIAG_Check
REAL :: sumt, sumv, sumd, value
INTEGER :: jl, jk, js
! compare new minus old with tendency
835: sumv = 0.0
DO js=1,nsp
DO jk=1,2
DO jl=1,nlev
value = svo(jl,jk,js)-pdovor(jl,jk,js)-SUM(pdvor(jl,jk,js,1:ndvor))
840: sumv = sumv + value*value
END DO
END DO
END DO
sumv = SQRT(sumv/(n2sp*nlev))
845:
sumd = 0.0
DO js=1,nsp
DO jk=1,2
DO jl=1,nlev
850: value = sd(jl,jk,js)-pdodiv(jl,jk,js)-SUM(pddiv(jl,jk,js,1:nddiv))
sumd = sumd + value*value
END DO
END DO
END DO
855: sumd = SQRT(sumd/(n2sp*nlev))
sumt = 0.0
DO js=1,nsp
DO jk=1,2
860: DO jl=1,nlev
value = stp(jl,jk,js)-pdotep(jl,jk,js)-SUM(pdtem(jl,jk,js,1:ndtem))
sumt = sumt + value*value
END DO
END DO
865: END DO
sumt = SQRT(sumt/(n2sp*nlev))
WRITE(nout,*) ' diagnose error VOR ',sumv,' DIV ',sumd,' TEM ',sumt
870: RETURN
END SUBROUTINE DIAG_Check
END MODULE mo_diag_tendency
Info Section
uses: mo_control, mo_doctor, mo_exception, mo_fft, mo_grib
mo_legendre, mo_memory_sp, mo_post, mo_start_dataset, mo_truncation
mo_year
calls: codegb5, diag_check, diag_init, diag_rerun, diagwritesp
fft991cy, finish, gribex, pbwrite, set_output_time
util_reshape
back to top
ECHAM 4 vf90 (C) 1998 Max-Planck-Institut für Meteorologie, Hamburg
Wed Nov 24 01:25:21 CST 1999
HTML derived from FORTRAN source by f2html.pl v0.3 (C) 1997,98 Beroud Jean-Marc.