si2.f90

      !+ 2nd part of the semi-implicit scheme (done in fourier space).
      !+ $Id: si2.f90,v 1.13 1999/07/22 13:15:04 m214030 Exp $
      !OCL NOALIAS
      
   5: SUBROUTINE si2
      
        ! Description:
        !
        ! 2nd part of the semi-implicit scheme (done in fourier space).
  10:   !
        ! Method:
        !
        ! This subroutine computes the contribution in
        ! *Fourier space to the semi-implicit scheme (mainly for
  15:   ! the vorticity and humidity equations).
        !
        ! *si2* is called from *fcc1*.
        !
        ! Reference:
  20:   ! 1-appendix *b1:organisation of the spectral model.
        !
        ! Authors:
        !
        ! M. Jarraud, ECMWF, January 1982, original source
  25:   ! L. Kornblueh, MPI, May 1998, f90 rewrite
        ! U. Schulzweida, MPI, May 1998, f90 rewrite
        ! I. Kirchner, MPI, October 1998, tendency diagnostics
        ! I. Kirchner, MPI, December 1998, vorticity semi-implicit diagnosticsterm modified
        ! 
  30:   ! for more details see file AUTHORS
        !
      
        USE mo_control,       ONLY: lptime, ltdiag, twodt
        USE mo_gaussgrid,     ONLY: racst, twomu
  35:   USE mo_truncation,    ONLY: am
        USE mo_semi_impl,     ONLY: betazq
        USE mo_constants,     ONLY: a
        USE mo_start_dataset, ONLY: nstart, nstep
        USE mo_diag_tendency, ONLY: pdiga, pdigb, pdigs
  40:   USE mo_decomposition, ONLY: dc=>local_decomposition
        USE mo_buffer_fft,    ONLY: fdm1, fvol, fvom, fu0, fdu0, fvom1
      
        IMPLICIT NONE
      
  45:   !  Local scalars: 
        REAL :: z1, z2, z3, zbmi, zbmr, zcmi, zcmr, zdl, zdt, zdtma2, zdtmda2, &
      &      zdu0, zrcst, zu0, zvoli, zvolr, zvomi, zvomr, zdigsr, zdigsi, &
      &      zdigsmr, zdigsmi
        INTEGER :: irow, jlev, jm, jrow, krow, jglat
  50:   INTEGER :: nflev, nmp1, ngl, nflat
      
        !  Executable statements 
      
        ! loop bounds on this PE
  55: 
        nflev = dc% nflev  ! local (Fourier space) : number of levels 
        nflat = dc% nflat  ! local (Fourier space) : number of latitudes  
        nmp1  = dc% nm+1   ! global                : number of coefficients m
        ngl   = dc% nlat   ! global                : number of latitudes
  60: 
      !-- 1. Locate and allocate storage
      
      !-- 1.1 Fourier components
      
  65:   ! latitude loop indices
      
        DO jrow = 1, nflat
      ! jrow  = nrow(2)                        ! local  continuous north to south
        jglat = dc% glat(jrow)                 ! global continuous north -> south
  70:   irow  = MIN(2*jglat-1,2*(ngl+1-jglat)) ! global ping pong index
        krow  = (irow+1)/2                     ! half latitude index
      
      !  ALLOCATE (dl(nlp2,nflev))
      
  75: !-- 2. Skip over *si2* during initialisation iterations
      !      or compute temporary quantities.
      
      !-- 2.1 Compute temporary quantities
      
  80:   zdt = twodt*.5
        IF (nstep==nstart) zdt = zdt*.5
        z1 = betazq*zdt*racst(krow)
        z2 = z1*a
        z3 = twomu(irow)*racst(krow)
  85:   zrcst = racst(krow)*a
      
      !-- 3. Semi implicit modifications
      
      !-- 3.1 Initiate scans
  90: 
        DO jlev = 1, nflev
          zu0  = z1*fu0(jlev,jrow)
          zdu0 = z2*(fdu0(jlev,jrow)+z3*fu0(jlev,jrow))
      
  95: !DIR$ IVDEP
      !OCL NOVREC
          DO jm = 1, nmp1
            zdtma2 = am(jm)*zu0
            zdtmda2 = am(jm)*zdu0
 100: 
            zdl = am(jm)*zrcst
      
            zbmr = 1./(1.+zdtma2**2)
            zbmi = -zdtma2*zbmr
 105: 
            zcmr = -2.*zdtmda2*    zdtma2    *zbmr**2
            zcmi =    -zdtmda2*(1.-zdtma2**2)*zbmr**2
      
      !-- 3.2 Divergence equation
 110: 
            fdm1(2*jm-1,jlev,jrow) =fdm1(2*jm-1,jlev,jrow)+zdl*fvom(2*jm,  jlev,jrow)
            fdm1(2*jm,  jlev,jrow) =fdm1(2*jm,  jlev,jrow)-zdl*fvom(2*jm-1,jlev,jrow)
      
      !-- 3.3 Vorticity equation
 115:       zvolr = zbmr*(fvom1(2*jm-1,jlev,jrow) - am(jm)*fvol(2*jm,  jlev,jrow))  &
      &           - zbmi*(fvom1(2*jm,  jlev,jrow) + am(jm)*fvol(2*jm-1,jlev,jrow))  &
      &           - zcmr* fvom (2*jm-1,jlev,jrow) &
      &           + zcmi* fvom (2*jm,  jlev,jrow)
      
 120:       zvoli = zbmi* (fvom1(2*jm-1,jlev,jrow) - am(jm)*fvol(2*jm,  jlev,jrow)) &
      &           + zbmr* (fvom1(2*jm,  jlev,jrow) + am(jm)*fvol(2*jm-1,jlev,jrow)) &
      &           - zcmi* fvom (2*jm-1,jlev,jrow) &
      &           - zcmr* fvom (2*jm,  jlev,jrow)
      
 125:       zvomr = zbmr*fvom(2*jm-1,jlev,jrow) - zbmi*fvom(2*jm,jlev,jrow)
            zvomi = zbmi*fvom(2*jm-1,jlev,jrow) + zbmr*fvom(2*jm,jlev,jrow)
      
            IF (ltdiag) THEN
               ! store M-Term before adjustment
 130:          zdigsmr = fvom(2*jm-1,jlev,jrow)
               zdigsmi = fvom(2*jm  ,jlev,jrow)
            ENDIF
      
            fvol(2*jm-1,jlev,jrow) = zvolr
 135:       fvol(2*jm,  jlev,jrow) = zvoli
            fvom(2*jm-1,jlev,jrow) = zvomr
            fvom(2*jm,  jlev,jrow) = zvomi
      
            IF (ltdiag) THEN
 140:          ! adjustment of semi-implicit part of vorticity, L-Term
               ! compose explicit part
               ! compose explicit part, accounted in SI1
               zdigsr = - am(jm)*pdigs(2*jm  ,jlev,1,irow)
               zdigsi =   am(jm)*pdigs(2*jm-1,jlev,1,irow)
 145:          ! add implicit part
      !         pdigs(2*jm-1,jlev,1,irow) = zdigsr + zdtma2*fvol(2*jm,  jlev,jrow)
      !         pdigs(2*jm,  jlev,1,irow) = zdigsi - zdtma2*fvol(2*jm-1,jlev,jrow)
               pdigs(2*jm-1,jlev,1,irow) = zdigsr + zu0*fvol(2*jm,  jlev,jrow)
               pdigs(2*jm,  jlev,1,irow) = zdigsi - zu0*fvol(2*jm-1,jlev,jrow)
 150:          ! adjustment of semi-implicit part of vorticity, M-Term
               pdigs(2*jm-1,jlev,4,irow) =   zdtma2*(zdigsmi-fvom(2*jm,  jlev,jrow))
               pdigs(2*jm,  jlev,4,irow) = - zdtma2*(zdigsmr-fvom(2*jm-1,jlev,jrow))
            ENDIF
      
 155:     END DO
      
          IF (ltdiag.AND.lptime) THEN
             ! zonal derivatives with factor m/(1-mu**2)
             DO jm = 1,nmp1
 160:           zdl = am(jm)*zrcst
                pdigb(2*jm-1,jlev,:,irow) = - zdl * pdiga(2*jm,  jlev,1:10,irow)
                pdigb(2*jm,  jlev,:,irow) =   zdl * pdiga(2*jm-1,jlev,1:10,irow)
             ENDDO
          ENDIF
 165:  
        END DO
      END DO
      END SUBROUTINE si2


Info Section
uses: mo_buffer_fft, mo_constants, mo_control, mo_decomposition, mo_diag_tendency mo_gaussgrid, mo_semi_impl, mo_start_dataset, mo_truncation
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.