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.