Changeset 2642 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Mar 15, 2022, 9:06:24 AM (3 years ago)
Author:
emillour
Message:

Mars GCM:
Switching orographic GW routines to F90 and adding comments.
JL

Location:
trunk/LMDZ.MARS
Files:
1 edited
6 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2640 r2642  
    36543654Update for non-orographic Gravity Waves; some key parameters can now be
    36553655set via callphys.def
     3656
     3657== 15/03/2022 == JL
     3658Switching orographic GW routines to F90 and adding comments.

  • trunk/LMDZ.MARS/libf/phymars/calldrag_noro_mod.F90

    r2641 r2642  
    1       MODULE calldrag_noro_mod
     1MODULE calldrag_noro_mod
     2     
     3IMPLICIT NONE
     4     
     5CONTAINS
     6     
     7      SUBROUTINE calldrag_noro(ngrid,nlayer,ptimestep,pplay,pplev,pt,pu,pv, &
     8                  !output
     9                  pdtgw,pdugw,pdvgw)
     10                 
     11!=====================================================================================================
     12! MODULE designed to call SUBROUTINE drag_noro_mod Interface for sub-grid scale orographic scheme
     13! The purpose of this subroutine is:
     14! 1) Make some initial calculation at first call
     15! 2) Split the calculation in several sub-grid ("sub-domain") to save memory and be able run on
     16!   a workstation at high resolution.The sub-grid size is defined in dimradmars_mod.
     17! 3) Transfer the increment output of drag_noro_mod into tendencies
     18! Christophe Hourdin + Francois Forget
     19!
     20! VERSION Update: This version uses the variable's splitting, which can be usefull when performing
     21! very high resolution simulation like LES. !
     22! J.-B. Madeleine 10W12
     23!
     24! UPDATE   Rewirten into .F90     J.Liu   3/3/2022      Translate the code into .F90. The name of the
     25!                                                       variables are uniformed. Comments are made.
     26!                                                       Unused variables are deleted.
     27!=======================================================================================================
     28
     29      use surfdat_h, only: zstd, zsig, zgam, zthe   !sub-grid scale standard devitation,slope,anisotropy, and priciple axes angle
     30                                                    ! which will be coded as pvar, psig,pgam,pthe in the called subroutines.
     31      use dimradmars_mod, only: ndomainsz
     32      use drag_noro_mod, only: drag_noro
     33      USE ioipsl_getin_p_mod, ONLY : getin_p
    234     
    335      IMPLICIT NONE
    4      
    5       CONTAINS
    6      
    7       SUBROUTINE calldrag_noro(ngrid,nlayer,ptimestep,
    8      &                 pplay,pplev,pt,pu,pv,pdtgw,pdugw,pdvgw)
    936
     37      ! 0.1  Declarations :
    1038
     39      ! 0.1.0 Input
     40      INTEGER, intent(in):: ngrid                 ! Number of atmospheric columns
     41      INTEGER, intent(in):: nlayer                ! Number of atmospheric layers
     42      REAL, intent(in):: ptimestep                ! Time step of the Physics(s)
     43      REAL, intent(in):: pplev(ngrid,nlayer+1)    ! Pressure at 1/2 levels(Pa)
     44      REAL, intent(in):: pplay(ngrid,nlayer)      ! Pressure at full levels(Pa)
     45      REAL, intent(in):: pt(ngrid,nlayer)         ! Temp at full levels(K)
     46      REAL, intent(in):: pu(ngrid,nlayer)         ! Zonal wind at full levels(m/s)
     47      REAL, intent(in):: pv(ngrid,nlayer)         ! Meridional wind at full levels(m/s)
     48     
     49      ! 0.1.1 Output   
     50      REAL, intent(out):: pdtgw(ngrid,nlayer)     ! Tendency on temperature (K/s) due to orographic gravity waves
     51      REAL, intent(out):: pdugw(ngrid,nlayer)     ! Tendency on zonal wind (m/s/s) due to orographic gravity waves
     52      REAL, intent(out):: pdvgw(ngrid,nlayer)     ! Tendency on meridional wind (m/s/s) due to orographic gravity waves
    1153
    12        use surfdat_h, only: zstd, zsig, zgam, zthe
    13        use dimradmars_mod, only: ndomainsz
    14        use drag_noro_mod, only: drag_noro
    15        IMPLICIT NONE
    16 c=======================================================================
    17 c   subject:
    18 c   --------
    19 c   Subroutine designed to call SUBROUTINE drag_noro
    20 c   Interface for sub-grid scale orographic scheme
    21 c   The purpose of this subroutine is
    22 c      1) Make some initial calculation at first call
    23 c      2) Split the calculation in several sub-grid
    24 c        ("sub-domain") to save memory and
    25 c        be able run on a workstation at high resolution
    26 c        The sub-grid size is defined in dimradmars_mod.
    27 c
    28 c   author:   
    29 c   ------
    30 c           Christophe Hourdin/ Francois Forget
    31 c
    32 c   changes:
    33 c   -------
    34 c   > J.-B. Madeleine 10W12
    35 c   This version uses the variable's splitting, which can be usefull
    36 c     when performing very high resolution simulation like LES.
    37 c
    38 c   input:
    39 c   -----
    40 c   ngrid                 number of gridpoint of horizontal grid
    41 c   nlayer                Number of layer
    42 c   ptimestep             Physical timestep (s)
    43 c   pplay(ngrid,nlayer)    pressure (Pa) in the middle of each layer
    44 c   pplev(ngrid,nlayer+1)  pressure (Pa) at boundaries of each layer
    45 c   pt(ngrid,nlayer)       atmospheric temperature  (K)
    46 c   pu(ngrid,nlayer)       zonal wind (m s-1)
    47 c   pv(ngrid,nlayer)       meridional wind (m s-1)
    48 c
    49 c   output:
    50 c   -------
    51 c   pdtgw(ngrid,nlayer)    Temperature trend (K.s-1)
    52 c   pdugw(ngrid,nlayer)    zonal wind trend  (m.s-2)
    53 c   pdvgw(ngrid,nlayer)    meridional wind trend  (m.s-2)
    54 c
    55 c
    56 c
    57 c
    58 c
    59 c=======================================================================
    60 c
    61 c    0.  Declarations :
    62 c    ------------------
    63 c
    64 
    65 c-----------------------------------------------------------------------
    66 c    Input/Output
    67 c    ------------
    68       INTEGER ngrid,nlayer 
    69 
    70       real ptimestep
    71 
    72       REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
    73       REAL pt(ngrid,nlayer), pu(ngrid,nlayer),pv(ngrid,nlayer)
    74       REAL pdtgw(ngrid,nlayer), pdugw(ngrid,nlayer),pdvgw(ngrid,nlayer)
    75 
    76 
    77 c
    78 c    Local variables :
    79 c    -----------------
    80 
    81       REAL sigtest(nlayer+1)
    82       INTEGER igwd,igwdim,itest(ngrid)
    83 
    84       INTEGER :: ndomain
    85 !      parameter (ndomain = (ngrid-1) / ndomainsz + 1)
    86 
     54      !0.1.2 Local variables :
     55      REAL sigtest(nlayer+1)      ! sigma coordinate at 1/2 level
     56      INTEGER kgwd                ! Number of points at which the scheme is called
     57      INTEGER kgwdim              ! kgwdIM=MAX(1,kgwd)
     58      INTEGER ktest(ngrid)        ! Map of calling points
     59      INTEGER kdx(ndomainsz)      ! Points at which to call the scheme
     60      INTEGER ndomain             ! Parameter (ndomain = (ngrid-1) / ndomainsz + 1)
    8761      INTEGER l,ig
    8862      INTEGER jd,ig0,nd
    89 
    90       REAL zulow(ngrid),zvlow(ngrid)
    91       REAL zustr(ngrid),zvstr(ngrid)
    92 
    93       REAL zplev(ndomainsz,nlayer+1)
    94       REAL zplay(ndomainsz,nlayer)
    95       REAL zt(ndomainsz,nlayer)
    96       REAL zu(ndomainsz,nlayer)
    97       REAL zv(ndomainsz,nlayer)
    98       INTEGER zidx(ndomainsz)
    99       REAL zzdtgw(ndomainsz,nlayer)
    100       REAL zzdugw(ndomainsz,nlayer)
    101       REAL zzdvgw(ndomainsz,nlayer)
    102 
    103       logical ll
    104 
    105 
    106 c   local saved variables
    107 c   ---------------------
    108 
    109       LOGICAL firstcall
    110 
     63      REAL,save:: zstdthread      ! Detecting points concerned by the scheme by height
     64!$OMP THREADPRIVATE(zstdthread)     
     65      REAL zulow(ngrid)              ! Low level zonal wind
     66      REAL zvlow(ngrid)              ! Low level Meridional wind
     67      REAL zustr(ngrid)              ! Low level zonal stress
     68      REAL zvstr(ngrid)              ! Low level meridional stress
     69      REAL zplev(ndomainsz,nlayer+1) ! pplev in sub domain
     70      REAL zplay(ndomainsz,nlayer)   ! pplay in sub domain
     71      REAL zt(ndomainsz,nlayer)      ! pt in sub domain
     72      REAL zu(ndomainsz,nlayer)      ! pu in sub domain
     73      REAL zv(ndomainsz,nlayer)      ! pv in sub domain
     74      REAL d_t(ndomainsz,nlayer)     ! Temperature increment(K) from DRAG_NORO
     75      REAL d_u(ndomainsz,nlayer)     ! Zonal wind increment(K=m/s) from DRAG_NORO
     76      REAL d_v(ndomainsz,nlayer)     ! Meridional wind increment(K=m/s) from DRAG_NORO
     77      logical,save:: ll=.false.
     78!$OMP THREADPRIVATE(ll)
     79      LOGICAL,save:: firstcall=.true.
    11180!$OMP THREADPRIVATE(firstcall)
    11281
    113       DATA firstcall/.true./
    114       SAVE firstcall
     82!-----------------------------------------------------------------------------------------------------------------------
     83!  1. INITIALISATIONS
     84!-----------------------------------------------------------------------------------------------------------------------
     85      IF (firstcall) THEN     
     86        write(*,*) "calldrag_noro: orographic GW scheme is active!"
     87        zstdthread = 50.0  ! Mars' value !!
     88        call getin_p("oro_gwd_zstdthread", zstdthread)
     89        write(*,*) "oro_gwd_zstdthread: zstdthread=", zstdthread
     90       
     91        do l=1,nlayer+1
     92           sigtest(l)=pplev(1,l)/pplev(1,1)
     93        enddo
     94         
     95        call sugwd(nlayer,sigtest)
    11596
    116 
    117 c----------------------------------------------------------------------
    118 
    119 c     Initialisation
    120 c     --------------
    121 
    122       IF (firstcall) THEN
    123 
    124          do l=1,nlayer+1
    125            sigtest(l)=pplev(1,l)/pplev(1,1)
    126          enddo
    127          call sugwd(nlayer,sigtest)
    128 
    129          if (ngrid .EQ. 1) then
    130            if (ndomainsz .NE. 1) then
     97        if (ngrid .EQ. 1) then
     98           IF (ndomainsz .NE. 1) THEN
    13199             print*
    132100             print*,'ATTENTION !!!'
    … …  
    135103             print*
    136104             call exit(1)
    137            endif
     105           ENDIF
     106        endif               
     107        firstcall=.false.
     108      END IF !firstcall
     109     
     110      ! Initialization the tendency of this scheme into zero
     111!      pdtgw(ngrid,nlayer)=0.0     ! Tendency on temperature (K/s) due to orographic gravity waves
     112!      pdugw(ngrid,nlayer)=0.0     ! Tendency on zonal wind (m/s/s) due to orographic gravity waves
     113!      pdvgw(ngrid,nlayer)=0.0     ! Tendency on meridional wind (m/s/s) due to orographic gravity waves
     114     
     115      ! AS: moved out of firstcall to allow nesting+evoluting horiz domain
     116      ndomain = (ngrid-1) / ndomainsz + 1
     117!-----------------------------------------------------------------------------------------------------------------------     
     118! 2. Starting loop on sub-domain
     119!-----------------------------------------------------------------------------------------------------------------------
     120      ! start the loop from the first sub-domain to the ndomain, in which the orographic gravity waves induced
     121      ! wind/temperature tendencies are calculated
     122      DO jd=1,ndomain
     123         ig0=(jd-1)*ndomainsz
     124         if (jd.eq.ndomain) then
     125            nd=ngrid-ig0
     126         else
     127            nd=ndomainsz
    138128         endif
     129                     
     130         ! 2.1 Detecting points concerned by the scheme
     131         kgwd=0
     132         DO ig=ig0+1,ig0+nd
     133            ktest(ig)=0
     134            ! only the points have a zstd greater than the thread value(default is 50) are concerned.
     135            ll=zstd(ig).gt.zstdthread
     136            IF(ll) then
     137               ! The map of calling is set 1 (which means the orographic map in this points is called?)
     138               ktest(ig)=1
     139               ! The numbers of the calling points added
     140               kgwd=kgwd+1
     141               ! The position of the calling points are picked up
     142               kdx(kgwd)=ig - ig0
     143            ENDIF
     144         ENDDO
     145         kgwdIM=MAX(1,kgwd)
    139146
    140          firstcall=.false.
    141       END IF
     147         ! 2.2 Spliting input variable in sub-domain input variables in each level
     148         DO l=1,nlayer+1
     149            do ig = 1,nd
     150               zplev(ig,l) = pplev(ig0+ig,l)
     151            enddo
     152         end DO
     153         DO l=1,nlayer
     154            do ig = 1,nd
     155               zplay(ig,l) = pplay(ig0+ig,l)
     156               zt(ig,l) = pt(ig0+ig,l)
     157               zu(ig,l) = pu(ig0+ig,l)
     158               zv(ig,l) = pv(ig0+ig,l)
     159            enddo
     160         end DO
    142161
    143       !! AS: moved out of firstcall to allow nesting+evoluting horiz domain
    144       ndomain = (ngrid-1) / ndomainsz + 1
     162         ! 2.3 Calling gravity wave and subgrid scale topo parameterization(all physics in this subroutine) to get
     163         ! zonal wind, meridional wind, and temperature increament
     164         call drag_noro (nd,nlayer,ptimestep,zplay,zplev,zstd(ig0+1),zsig(ig0+1),zgam(ig0+1),zthe(ig0+1),&
     165              kgwd,kgwdim,kdx,ktest(ig0+1),zt, zu, zv, &
     166              ! output
     167              zulow(ig0+1),zvlow(ig0+1),zustr(ig0+1),zvstr(ig0+1),d_t,d_u,d_v)
    145168
    146 c     Starting loop on sub-domain
    147 c     ----------------------------
     169         ! 2.4  Un-spliting output variable from sub-domain input variables (and devide by ptimestep -> true tendencies)
     170         DO l=1,nlayer
     171            do ig = 1,nd
     172                   pdtgw(ig0+ig,l) = d_t(ig,l)/ptimestep
     173                   pdugw(ig0+ig,l) = d_u(ig,l)/ptimestep
     174                   pdvgw(ig0+ig,l) = d_v(ig,l)/ptimestep
     175            enddo
     176         end DO
    148177
    149       DO jd=1,ndomain
    150         ig0=(jd-1)*ndomainsz
    151         if (jd.eq.ndomain) then
    152           nd=ngrid-ig0
    153         else
    154           nd=ndomainsz
    155         endif
     178      end DO ! (jd=1, ndomain)     
     179!     2. ending loop on sub-domain
    156180
    157 c       Detecting points concerned by the scheme
    158 c       ----------------------------------------
     181!******************************************************************************
    159182
    160         igwd=0
    161         DO ig=ig0+1,ig0+nd
    162           itest(ig)=0
    163           ll=zstd(ig).gt.50.0
    164           IF(ll) then
    165             itest(ig)=1
    166             igwd=igwd+1
    167             zidx(igwd)=ig - ig0
    168           ENDIF
    169         ENDDO
    170         IGWDIM=MAX(1,IGWD)
     183      END SUBROUTINE calldrag_noro  ! end of the subroutine   
     184     
     185END MODULE calldrag_noro_mod  ! end of the module
    171186
    172 c       Spliting input variable in sub-domain input variables
    173 c       ---------------------------------------------------
    174 
    175         do l=1,nlayer+1
    176           do ig = 1,nd
    177            zplev(ig,l) = pplev(ig0+ig,l)
    178           enddo
    179         enddo
    180 
    181         do l=1,nlayer
    182           do ig = 1,nd
    183             zplay(ig,l) = pplay(ig0+ig,l)
    184             zt(ig,l) = pt(ig0+ig,l)
    185             zu(ig,l) = pu(ig0+ig,l)
    186             zv(ig,l) = pv(ig0+ig,l)
    187           enddo
    188         enddo
    189 
    190 c       Calling gravity wave and subgrid scale topo parameterization
    191 c       -------------------------------------------------------------
    192 
    193         call drag_noro (nd,nlayer,ptimestep,zplay,zplev,
    194      e        zstd(ig0+1),zsig(ig0+1),zgam(ig0+1),zthe(ig0+1),
    195      e        igwd,igwdim,zidx,itest(ig0+1),
    196      e        zt, zu, zv,
    197      s        zulow(ig0+1),zvlow(ig0+1),zustr(ig0+1),zvstr(ig0+1),
    198      s        zzdtgw,zzdugw,zzdvgw)
    199 
    200 c       Un-spliting output variable from sub-domain input variables
    201 c       ------------------------------------------------------------
    202 c       (and devide by ptimestep -> true tendancies)
    203 
    204         do l=1,nlayer
    205          do ig = 1,nd
    206           pdtgw(ig0+ig,l) = zzdtgw(ig,l)/ptimestep
    207           pdugw(ig0+ig,l) = zzdugw(ig,l)/ptimestep
    208           pdvgw(ig0+ig,l) = zzdvgw(ig,l)/ptimestep
    209          enddo
    210         enddo
    211 
    212       ENDDO         !   (boucle jd=1, ndomain)
    213 
    214       END SUBROUTINE calldrag_noro
    215      
    216       END MODULE calldrag_noro_mod
    217 

  • trunk/LMDZ.MARS/libf/phymars/drag_noro_mod.F90

    r2641 r2642  
    1       MODULE drag_noro_mod
     1MODULE drag_noro_mod
    22     
    3       IMPLICIT NONE
     3IMPLICIT NONE
    44     
    5       CONTAINS
     5CONTAINS
    66     
    7       SUBROUTINE drag_noro (klon,klev,dtime,pplay,pplev,
    8      e                   pvar, psig, pgam, pthe,
    9      e                   kgwd,kgwdim,kdx,ktest,
    10      e                   t, u, v,
    11      s                   pulow, pvlow, pustr, pvstr,
    12      s                   d_t, d_u, d_v)
    13 C**** *DRAG_NORO* INTERFACE FOR SUB-GRID SCALE OROGRAPHIC SCHEME
    14 C
    15 C     PURPOSE.
    16 C     --------
    17 C           ZEROS TENDENCIES, COMPUTES GEOPOTENTIAL HEIGHT AND UPDATES THE
    18 C           TENDENCIES AFTER THE SCHEME HAS BEEN CALLED.
    19 C
    20 C     EXPLICIT ARGUMENTS :
    21 C     --------------------
    22 C
    23 C     INPUT :
    24 C
    25 C     NLON               : NUMBER OF HORIZONTAL GRID POINTS
    26 C     NLEV               : NUMBER OF LEVELS
    27 C     DTIME              : LENGTH OF TIME STEP
    28 C     PPLAY(NLON,NLEV+1) : PRESSURE AT MIDDLE LEVELS
    29 C     PPLEV(NLON,NLEV)   : PRESSURE ON MODEL LEVELS
    30 C     PVAR(NLON)         : SUB-GRID SCALE STANDARD DEVIATION
    31 C     PSIG(NLON)         : SUB-GRID SCALE SLOPE
    32 C     PGAM(NLON)         : SUB-GRID SCALE ANISOTROPY
    33 C     PTHE(NLON)         : SUB-GRID SCALE PRINCIPAL AXES ANGLE
    34 C     KGWD               : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
    35 C     KGWDIM             : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
    36 C     KDX(NLON)          : POINTS AT WHICH TO CALL THE SCHEME
    37 C     KTEST(NLON)        : MAP OF CALLING POINTS
    38 C     T(NLON,NLEV)       : TEMPERATURE
    39 C     U(NLON,NLEV)       : ZONAL WIND
    40 C     V(NLON,NLEV)       : MERIDIONAL WIND
    41 C
    42 C     OUTPUT :
    43 C
    44 C     PULOW(NLON)        : LOW LEVEL ZONAL WIND
    45 C     PVLOW(NLON)        : LOW LEVEL MERIDIONAL WIND
    46 C     PUSTR(NLON)        : LOW LEVEL ZONAL STRESS
    47 C     PVSTR(NLON)        : LOW LEVEL MERIDIONAL STRESS
    48 C     D_T(NLON,NLEV)     : TEMPERATURE TENDENCY
    49 C     D_U(NLON,NLEV)     : ZONAL WIND TENDENCY
    50 C     D_V(NLON,NLEV)     : MERIDIONAL WIND TENDENCY
    51 C
    52 C     IMPLICIT ARGUMENTS :
    53 C     --------------------
    54 C
    55 C     comcstfi.h
    56 C
    57 c
    58       use dimradmars_mod, only:  ndlo2
     7      SUBROUTINE drag_noro (ngrid,nlayer,ptimestep,pplay,pplev,pvar, psig, pgam, pthe, &
     8                 kgwd,kgwdim,kdx,ktest,t, u, v, &
     9                 ! output
     10                 pulow, pvlow, pustr, pvstr, d_t, d_u, d_v)
     11
     12!--------------------------------------------------------------------------------
     13! MODULE contains DRAG_NORO SUBROUNTINE for sub-grid scale orographic sheme.
     14!  1) Initialization the variables
     15!  2) Overturn the levels and computes geopotential height
     16!  3) Call the and ORODRAG subroutine to updates the tendencies
     17!  4) Overturn back to the normal level of the tendencies and transfer it into
     18!     increments and sum the oro-gw stress.
     19! the scheme has been called.         
     20! Z.X. Li + F.Lott (LMD/CNRS/ENS)  1995-02-01
     21!
     22! UPDATE: J.Liu   03/03/2022      Translate the code into .F90. The name of the
     23!                                 variables are uniformed. Comments are made.
     24!                                 Unused variables are deleted.
     25!-------------------------------------------------------------------------------
     26 
     27      use dimradmars_mod, only:  ndomainsz
    5928      USE orodrag_mod, ONLY: orodrag
    6029      USE comcstfi_h, ONLY: g, r
     30     
    6131      IMPLICIT none
    62 c======================================================================
    63 c Auteur(s): Z.X. Li F.Lott (LMD/CNRS) date: 19950201
    64 c Objet: Frottement de la montagne Interface
    65 c======================================================================
    66 c Arguments:
    67 c dtime---input-R- pas d'integration (s)
    68 c s-------input-R-la valeur "s" pour chaque couche
    69 c pplay--input-R- pression au milieu des couches en Pa
    70 c pplev--input-R-pression au bords des couches en Pa
    71 c t-------input-R-temperature (K)
    72 c u-------input-R-vitesse horizontale (m/s)
    73 c v-------input-R-vitesse horizontale (m/s)
    74 c
    75 c d_t-----output-R-increment de la temperature t
    76 c d_u-----output-R-increment de la vitesse u
    77 c d_v-----output-R-increment de la vitesse v
    78 c======================================================================
    79 c
    80 c ARGUMENTS
    81 c
    82       REAL dtime
    83       INTEGER klon,klev
    84       real pplay(NDLO2,klev),pplev(NDLO2,klev+1)
    85       REAL pvar(NDLO2),psig(NDLO2),pgam(NDLO2),pthe(NDLO2)
    86       REAL pulow(NDLO2),pvlow(NDLO2),pustr(NDLO2),pvstr(NDLO2)
    87       REAL u(NDLO2,klev), v(NDLO2,klev),t(NDLO2,klev)
    88       REAL d_t(NDLO2,klev), d_u(NDLO2,klev), d_v(NDLO2,klev)
    89 c
    90       INTEGER i, k, kgwd, kgwdim, kdx(NDLO2), ktest(NDLO2)
    91 c
    92 c Variables locales:
    93 c
    94       REAL paprs(NDLO2,klev+1)
    95       REAL paprsf(NDLO2,klev)
    96       REAL zgeom(NDLO2,klev)
    97       REAL pdtdt(NDLO2,klev)
    98       REAL pdudt(NDLO2,klev), pdvdt(NDLO2,klev)
    99       REAL pt(NDLO2,klev), pu(NDLO2,klev)
    100       REAL pv(NDLO2,klev)
    101 c
    102 c initialiser les variables de sortie (pour securite)
    103 c
    104       DO i = 1,klon
     32     
     33      ! 0. DECLARATIONS:
     34     
     35      ! 0.1 Input:
     36      INTEGER, intent(in):: ngrid                       ! Number of atmospheric columns [only nd]
     37      INTEGER, intent(in):: nlayer                      ! Number of atmospheric layers
     38      REAL, intent(in):: ptimestep                      ! Time step of the Physics(s)
     39      real, intent(in):: pplay(ndomainsz,nlayer)        ! Pressure at full levels(Pa)
     40      real, intent(in):: pplev(ndomainsz,nlayer+1)      ! Pressure at 1/2 levels(Pa)
     41      REAL, intent(in):: pvar(ndomainsz)                ! sub-grid scale standard deviation
     42      REAL, intent(in):: psig(ndomainsz)                ! sub-grid scale standard slope
     43      REAL, intent(in):: pgam(ndomainsz)                ! sub-grid scale standard anisotropy
     44      REAL, intent(in):: pthe(ndomainsz)                ! sub-grid scale principal axes angle
     45      REAL, intent(in):: u(ndomainsz,nlayer)            ! Zonal wind at full levels(m/s)
     46      REAL, intent(in):: v(ndomainsz,nlayer)            ! Meridional wind at full levels(m/s)
     47      REAL, intent(in):: t(ndomainsz,nlayer)            ! Temperature at full levels(m/s)
     48      INTEGER, intent(in):: kgwd                        ! Number of points at which the scheme is called
     49      INTEGER, intent(in):: kgwdim                      ! kgwdim=MAX(1,kgwd)
     50      INTEGER, intent(in):: kdx(ndomainsz)              ! Points at which to call the scheme
     51      INTEGER, intent(in):: ktest(ndomainsz)            ! Map of calling points
     52     
     53      ! 0.2 Output:
     54      REAL, intent(out):: pulow(ndomainsz)              ! Low level zonal wind
     55      REAL, intent(out):: pvlow(ndomainsz)              ! Low level meridional wind
     56      REAL, intent(out):: pustr(ndomainsz)              ! Low level zonal stress     
     57      REAL, intent(out):: pvstr(ndomainsz)              ! Low level meridional stress 
     58      REAL, intent(out):: d_t(ndomainsz,nlayer)         ! Temperature increment(K) due to orographic gravity waves     
     59      REAL, intent(out):: d_u(ndomainsz,nlayer)         ! Zonal increment(m/s) due to orographic gravity waves
     60      REAL, intent(out):: d_v(ndomainsz,nlayer)         ! Meridional increment(m/s) due to orographic gravity waves
     61     
     62      ! 0.3 Variables locales:
     63      INTEGER i, k
     64      REAL inv_pplev(ndomainsz,nlayer+1)                ! Inversed (by inverse the pplev pressure) pressure at 1/2 levels
     65      REAL inv_pplay(ndomainsz,nlayer)                  ! Inversed (by inverse the pplay pressure) pressure at full levels
     66      REAL zgeom(ndomainsz,nlayer)                      ! Geopotetial height (Inversed ??)
     67      REAL inv_pt(ndomainsz,nlayer)                     ! Inversed (by inverse the t) temperature (K) at full levels
     68      REAL inv_pu(ndomainsz,nlayer)                     ! Inversed (by inverse the u) zonal wind (m/s) at full levels
     69      REAL inv_pv(ndomainsz,nlayer)                     ! Inversed (by inverse the v) meridional wind (m/s) at full levels
     70      REAL pdtdt(ndomainsz,nlayer)                      ! Temperature tendency outputs from main routine ORODRAG
     71      REAL pdudt(ndomainsz,nlayer)                      ! Zonal winds tendency outputs from main routine ORODRAG
     72      REAL pdvdt(ndomainsz,nlayer)                      ! Meridional winds tendency outputs from main routine ORODRAG
     73     
     74!-----------------------------------------------------------------------------------------------------------------------
     75!  1. INITIALISATIONS
     76!-----------------------------------------------------------------------------------------------------------------------
     77      ! 1.1 Initialize the input/output variables (for security purpose)
     78      DO i = 1,ngrid
    10579         pulow(i) = 0.0
    10680         pvlow(i) = 0.0
    … …  
    10882         pvstr(i) = 0.0
    10983      ENDDO
    110       DO k = 1, klev
    111       DO i = 1, klon
     84     
     85      DO k = 1, nlayer
     86        DO i = 1, ngrid
    11287         d_t(i,k) = 0.0
    11388         d_u(i,k) = 0.0
    … …  
    11691         pdvdt(i,k)=0.0
    11792         pdtdt(i,k)=0.0
     93        ENDDO
    11894      ENDDO
     95     
     96      ! 1.2 Prepare the ininv_put varibales (Attention: The order of vertical levels increases from top to bottom )
     97      ! Here the levels are overturned, that is, the surface becomes to the top and the top becomes to the surface
     98      ! and so forth
     99      DO k = 1, nlayer
     100        DO i = 1, ngrid
     101         inv_pt(i,k) = t(i,nlayer-k+1)
     102         inv_pu(i,k) = u(i,nlayer-k+1)
     103         inv_pv(i,k) = v(i,nlayer-k+1)
     104         inv_pplay(i,k) = pplay(i,nlayer-k+1)
     105         inv_pplev(i,k) = pplev(i,nlayer+1-k+1)
     106        ENDDO
     107      endDO
     108           
     109      DO i = 1, ngrid
     110         inv_pplev(i,nlayer+1) = pplev(i,1)
    119111      ENDDO
    120 c
    121 c preparer les variables d'entree (attention: l'ordre des niveaux
    122 c verticaux augmente du haut vers le bas)
    123 c
    124       DO k = 1, klev
    125       DO i = 1, klon
    126          pt(i,k) = t(i,klev-k+1)
    127          pu(i,k) = u(i,klev-k+1)
    128          pv(i,k) = v(i,klev-k+1)
    129          paprsf(i,k) = pplay(i,klev-k+1)
    130          paprs(i,k) = pplev(i,klev+1-k+1)
     112     
     113      !calculate g*dz by R*T*[LN(P(z)/P(z+dz))]
     114      DO i = 1, ngrid
     115         zgeom(i,nlayer) = r * inv_pt(i,nlayer)* LOG(inv_pplev(i,nlayer+1)/inv_pplay(i,nlayer))
    131116      ENDDO
     117     
     118      ! sum g*dz from surface to top to get geopotential height
     119      DO k = nlayer-1, 1, -1
     120        DO i = 1, ngrid
     121           zgeom(i,k) = zgeom(i,k+1) + r * (inv_pt(i,k)+inv_pt(i,k+1))/2.0 * LOG(inv_pplay(i,k+1)/inv_pplay(i,k))
     122        ENDDO
     123      endDO
     124     
     125!-----------------------------------------------------------------------------------------------------------------------
     126!  2. CALL the Main Rountine OORDRAG
     127!-----------------------------------------------------------------------------------------------------------------------
     128      ! 2.1 call the main rountine to get the tendencies   
     129      CALL ORODRAG(ngrid,nlayer,kgwd,kgwdim,kdx,ktest,ptimestep,inv_pplev,inv_pplay,zgeom,&
     130           inv_pt, inv_pu, inv_pv, pvar, psig, pgam, pthe,                                &
     131           ! output:
     132           pulow,pvlow,pdudt,pdvdt,pdtdt)
     133                 
     134      ! 2.2 Transfer tendency into increment by multiply the physical time steps
     135      !     maybe in the future here cancel multiply the ptimestep thus it will not divide ptimestep again in the main rountine
     136      !     calldrag_noro_mod.F90
     137      DO k = 1, nlayer
     138        DO i = 1, ngrid
     139             d_u(i,nlayer+1-k) = ptimestep*pdudt(i,k)
     140             d_v(i,nlayer+1-k) = ptimestep*pdvdt(i,k)
     141             d_t(i,nlayer+1-k) = ptimestep*pdtdt(i,k)
     142             pustr(i)          = pustr(i)+g*pdudt(i,k)*(inv_pplev(i,k+1)-inv_pplev(i,k))
     143             pvstr(i)          = pvstr(i)+g*pdvdt(i,k)*(inv_pplev(i,k+1)-inv_pplev(i,k))
     144        ENDDO
    132145      ENDDO
    133       DO i = 1, klon
    134          paprs(i,klev+1) = pplev(i,1)
    135       ENDDO
    136       DO i = 1, klon
    137          zgeom(i,klev) = r * pt(i,klev)
    138      .                  * LOG(paprs(i,klev+1)/paprsf(i,klev))
    139       ENDDO
    140       DO k = klev-1, 1, -1
    141       DO i = 1, klon
    142          zgeom(i,k) = zgeom(i,k+1) + r * (pt(i,k)+pt(i,k+1))/2.0
    143      .               * LOG(paprsf(i,k+1)/paprsf(i,k))
    144       ENDDO
    145       ENDDO
    146 c
    147 c appeler la routine principale
    148 c
    149 
    150       CALL ORODRAG(klon,klev,kgwd,kgwdim,kdx,ktest,
    151      .            dtime,
    152      .            paprs, paprsf, zgeom,
    153      .            pt, pu, pv, pvar, psig, pgam, pthe,
    154      .            pulow,pvlow,
    155      .            pdudt,pdvdt,pdtdt)
    156 C
    157       DO k = 1, klev
    158       DO i = 1, klon
    159          d_u(i,klev+1-k) = dtime*pdudt(i,k)
    160          d_v(i,klev+1-k) = dtime*pdvdt(i,k)
    161          d_t(i,klev+1-k) = dtime*pdtdt(i,k)
    162          pustr(i)        = pustr(i)
    163      .                    +g*pdudt(i,k)*(paprs(i,k+1)-paprs(i,k))
    164          pvstr(i)        = pvstr(i)
    165      .                    +g*pdvdt(i,k)*(paprs(i,k+1)-paprs(i,k))
    166       ENDDO
    167       ENDDO
    168 c
    169146
    170147      END SUBROUTINE drag_noro
    171148     
    172       END MODULE drag_noro_mod
     149END MODULE drag_noro_mod

  • trunk/LMDZ.MARS/libf/phymars/gwprofil_mod.F90

    r2641 r2642  
    1       MODULE gwprofil_mod
     1MODULE gwprofil_mod
    22     
    3       IMPLICIT NONE
     3IMPLICIT NONE
    44     
    5       CONTAINS
     5CONTAINS
    66     
    7       SUBROUTINE GWPROFIL
    8      *         ( klon, klev
    9      *         , kgwd ,kdx  , ktest
    10      *         , KKCRIT, KKCRITH, KCRIT ,  kkenvh, kknu,kknu2
    11      *         , PAPHM1, PRHO   , PSTAB , PTFR , PVPH , PRI , PTAU
    12      *         , ptauf ,pdmod   , pnu   , psig ,pgamma, pvar      )
     7      SUBROUTINE GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, &
     8                 IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev,  &
     9                 ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar,    &
     10                 !not used variables
     11                 !IKCRIT,ZTAUf,ZTFR,znu,pgam,
     12                 ! in-output
     13                 ZTAU      )
    1314
    14 C**** *GWPROFIL*
    15 C
    16 C     PURPOSE.
    17 C     --------
    18 C
    19 C**   INTERFACE.
    20 C     ----------
    21 C          FROM *GWDRAG*
    22 C
    23 C        EXPLICIT ARGUMENTS :
    24 C        --------------------
    25 C     ==== INPUTS ===
    26 C     ==== OUTPUTS ===
    27 C
    28 C        IMPLICIT ARGUMENTS :   NONE
    29 C        --------------------
    30 C
    31 C     METHOD:
    32 C     -------
    33 C     THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
    34 C     IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
    35 C     AND THE TOP OF THE BLOCKED LAYER (KKENVH).
    36 C     IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE
    37 C     BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR
    38 C     NONLINEAR GRAVITY WAVE BREAKING.
    39 C     ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
    40 C     LEVEL (KCRIT) OR WHEN IT BREAKS.
    41 C     
    42 C
    43 C
    44 C     EXTERNALS.
    45 C     ----------
    46 C
    47 C
    48 C     REFERENCE.
    49 C     ----------
    50 C
    51 C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
    52 C
    53 C     AUTHOR.
    54 C     -------
    55 C
    56 C     MODIFICATIONS.
    57 C     --------------
    58 C     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
    59 C-----------------------------------------------------------------------
    60       use dimradmars_mod, only: ndlo2
     15        !---------------------------------------------------------------------------
     16        !     THe stress profile for gravity waves is computed as follows:
     17        !     It is constant (no GWD) at the levels between the ground and
     18        !     the top of the  blocked layer (IKENVH). It decreases linearly
     19        !     with heights from the top of the blocked layer to 3*varor(IKNU)
     20        !     to simulates lee waves or nonlinear gravity wave breaking.
     21        !     Above it is constant, except when the wave encounters a critical
     22        !     level(ICRIT) or when it breaks.!
     23        !     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
     24        !     REFERENCE.
     25        !     SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
     26        !     MODIFICATIONS.
     27        !     Rewrited by J.Liu 03/03/2022
     28        !---------------------------------------------------------------------------
     29       
     30      use dimradmars_mod, only: ndomainsz 
    6131      implicit none
    62 C
     32#include "yoegwd.h"   ! why begin with #, needs ask Ehouarn Millour
     33 
     34      ! 0. DECLARATIONS:
     35      ! 0.1   ARGUMENTS
     36      integer, intent(in):: ngrid    ! Number of atmospheric columns
     37      integer, intent(in):: nlayer   ! Number of atmospheric layers
     38      INTEGER, intent(in) :: kgwd    ! Number of points at which the scheme is called     
     39!      INTEGER, INTENT(IN) :: IKCRIT(ndomainsz  ) ! Top of low level flow height (not used)
     40      INTEGER, INTENT(IN) :: IKCRITH(ndomainsz  ) ! dynamical mixing height for the breaking of gravity waves
     41      INTEGER, INTENT(IN) :: ICRIT(ndomainsz  )   ! Critical layer where orographic GW breaks
     42      INTEGER, INTENT(IN) :: kdx(ndomainsz  )     ! Points at which to call the scheme   
     43      INTEGER, INTENT(IN) :: ktest(ndomainsz  )   ! Map of calling points
     44      INTEGER, INTENT(IN) :: IKENVH(ndomainsz  )  ! Top of the  blocked layer
     45      INTEGER, INTENT(IN) :: IKNU(ndomainsz  )    ! 4*pvar layer
     46      INTEGER, INTENT(IN) :: IKNU2(ndomainsz  )   ! 3*pvar layer
    6347
    64 C
     48      REAL, INTENT(IN):: pplev(ndomainsz  ,nlayer+1)   ! Pressure at 1/2 levels(Pa)
     49      REAL, INTENT(IN):: BV(ndomainsz  ,nlayer+1)      ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL)
     50      REAL, INTENT(IN):: ZRHO  (ndomainsz  ,nlayer+1)  ! Atmospheric density at 1/2 level
     51      REAL, INTENT(IN):: ZVPH (ndomainsz  ,nlayer+1)   ! SUB-GRID SCALE STANDARD DEVIATION at 1/2 level
     52      REAL, INTENT(IN):: PRI   (ndomainsz  ,nlayer+1)  ! Mean flow richardson number at 1/2 level
     53!      REAL, INTENT(IN):: ZTFR (ndomainsz  )           ! not used
     54!      REAL, INTENT(IN):: ZTAUf (ndomainsz  ,nlayer+1) ! not used
     55!     
     56      REAL, INTENT(IN):: zdmod (ndomainsz  )
     57!      REAL, INTENT(IN):: znu (ndomainsz  )            ! not used
     58      REAL, INTENT(IN):: psig(ndomainsz  )             ! Sub-grid scale slope
     59!      REAL, INTENT(IN):: pgam(ndomainsz  )            ! Sub-grid scale anisotropy(not used)
     60      REAL, INTENT(IN):: pvar(ndomainsz  )             ! Sub-grid scale standard deviation
     61      REAL, INTENT(inout):: ZTAU(ndomainsz  ,nlayer+1) ! Gravity waves induced stress
    6562
    66       integer klon,klev,kidia,kfdia
    67 #include "yoegwd.h"
     63      ! 0.2   LOCAL ARRAYS
     64      real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt
     65      integer ji,jk,jl,ilevh
     66      REAL ZDZ2 (ndomainsz  ,nlayer) , ZNORM(ndomainsz  ) , zoro(ndomainsz  )
     67      REAL Z_TAU (ndomainsz  ,nlayer+1)
    6868
    69 C-----------------------------------------------------------------------
    70 C
    71 C*       0.1   ARGUMENTS
    72 C              ---------
    73 C
    74       integer kgwd
    75       INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2)
    76      *       ,kdx(NDLO2),ktest(NDLO2)
    77      *       ,kkenvh(NDLO2),kknu(NDLO2),kknu2(NDLO2)
    78 C
    79       REAL PAPHM1(NDLO2,klev+1), PSTAB(NDLO2,klev+1),
    80      *     PRHO  (NDLO2,klev+1), PVPH (NDLO2,klev+1),
    81      *     PRI   (NDLO2,klev+1), PTFR (NDLO2), PTAU(NDLO2,klev+1),
    82      *     ptauf (NDLO2,klev+1)
     69!-----------------------------------------------------------------------
     70!*         1.    Initialization
     71!-----------------------------------------------------------------------
     72!      kidia=1
     73!      kfdia=ngrid
     74 100  CONTINUE
     75      ! 1.1 Computional constants .
     76      ilevh=nlayer/3
     77
     78      DO ji=1,kgwd
     79        jl=kdx(ji)
     80        Zoro(JL)=psig(JL)*zdmod(JL)/4./max(pvar(jl),1.0)
     81        Z_TAU(JL,IKNU(JL)+1)=ZTAU(JL,IKNU(JL)+1)
     82        Z_TAU(JL,nlayer+1)=ZTAU(JL,nlayer+1)
     83      ENDDO
     84     
     85      ! all start here with this long loop
     86      DO JK=nlayer,2,-1   
     87        ! 4.1   Constant wave stress until top of the blocking layer
     88  410 CONTINUE
     89        DO ji=1,kgwd
     90           jl=kdx(ji)
     91           IF(JK.GE.IKNU2(JL)) THEN
     92           ZTAU(JL,JK)=Z_TAU(JL,nlayer+1)
     93           ENDIF
     94        ENDDO           
     95
     96      !*         4.15  Constant shear stress until the top of the low level flow layer
     97 415  CONTINUE
     98      ! 4.2    Wave displacement at next level.
     99  420 CONTINUE
     100
     101        DO ji=1,kgwd
     102           jl=kdx(ji)
     103           IF(JK.LT.IKNU2(JL)) THEN
     104                ZNORM(JL)=gkdrag*ZRHO(JL,JK)*SQRT(BV(JL,JK))*ZVPH(JL,JK)*zoro(jl)
     105                ZDZ2(JL,JK)=ZTAU(JL,JK+1)/max(ZNORM(JL),gssec)
     106           ENDIF
     107        ENDDO
     108
     109       ! 4.3    Wave Richardson number, new wave displacement and stress:
     110       ! Breaking evaluation and critical level C
     111        DO ji=1,kgwd
     112           jl=kdx(ji)
     113           IF(JK.LT.IKNU2(JL)) THEN
     114                IF((ZTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.ICRIT(JL))) THEN
     115                   ZTAU(JL,JK)=0.0
     116                ELSE
     117                   ZSQR=SQRT(PRI(JL,JK))
     118                   ZALFA=SQRT(BV(JL,JK)*ZDZ2(JL,JK))/ZVPH(JL,JK)
     119                   ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
     120                   IF(ZRIW.LT.GRCRIT) THEN
     121                     ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
     122                     ZB=1./GRCRIT+2./ZSQR
     123                     ZALPHA=0.5*(-ZB+SQRT(ZDEL))
     124                     ZDZ2N=(ZVPH(JL,JK)*ZALPHA)**2/BV(JL,JK)
     125                     ZTAU(JL,JK)=ZNORM(JL)*ZDZ2N
     126                   ELSE
     127                     ZTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
     128                   ENDIF
     129                   ZTAU(JL,JK)=MIN(ZTAU(JL,JK),ZTAU(JL,JK+1))
     130                ENDIF
     131           ENDIF
     132        ENDDO 
     133      end DO !JK=nlayer,2,-1
     134!-------------------------------------------------------------------------------------------
    83135     
    84       REAL pdmod (NDLO2) , pnu (NDLO2) , psig(NDLO2),
    85      *     pgamma(NDLO2) , pvar(NDLO2)
    86      
    87 C-----------------------------------------------------------------------
    88 C
    89 C*       0.2   LOCAL ARRAYS
    90 C              ------------
    91 C
    92 c   declarations pour 'implicit none"
    93       real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt
    94 
    95       integer ji,jk,jl,ilevh
    96       REAL ZDZ2 (NDLO2,klev) , ZNORM(NDLO2) , zoro(NDLO2)
    97       REAL ZTAU (NDLO2,klev+1)
    98 C
    99 C-----------------------------------------------------------------------
    100 C
    101 C*         1.    INITIALIZATION
    102 C                --------------
    103 
    104 
    105       kidia=1
    106       kfdia=klon
    107 
    108  100  CONTINUE
    109 C
    110 C
    111 C*    COMPUTATIONAL CONSTANTS.
    112 C     ------------- ----------
    113 C
    114       ilevh=KLEV/3
    115 C
    116       DO 400 ji=1,kgwd
    117       jl=kdx(ji)
    118       Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
    119       ZTAU(JL,KKNU(JL)+1)=PTAU(JL,KKNU(JL)+1)
    120       ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
    121   400 CONTINUE
    122 C
    123       DO 430 JK=KLEV,2,-1
    124 C
    125 C
    126 C*         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE
    127 C                 BLOCKING LAYER.
    128   410 CONTINUE
    129 C
    130       DO 411 ji=1,kgwd
    131       jl=kdx(ji)
    132            IF(JK.GE.KKNU2(JL)) THEN
    133            PTAU(JL,JK)=ZTAU(JL,KLEV+1)
    134            ENDIF
    135  411  CONTINUE             
    136 C
    137 C*         4.15   CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
    138 C                 LOW LEVEL FLOW LAYER.
    139  415  CONTINUE
    140 C       
    141 C
    142 C*         4.2    WAVE DISPLACEMENT AT NEXT LEVEL.
    143 C
    144   420 CONTINUE
    145 C
    146       DO 421 ji=1,kgwd
    147       jl=kdx(ji)
    148       IF(JK.LT.KKNU2(JL)) THEN
    149       ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
    150      *                                                    *zoro(jl)
    151       ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
    152       ENDIF
    153   421 CONTINUE
    154 C
    155 C*         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
    156 C*                AND STRESS:  BREAKING EVALUATION AND CRITICAL
    157 C                 LEVEL
    158 C
    159       DO 431 ji=1,kgwd
    160       jl=kdx(ji)
    161           IF(JK.LT.KKNU2(JL)) THEN
    162           IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
    163             PTAU(JL,JK)=0.0
    164           ELSE
    165                ZSQR=SQRT(PRI(JL,JK))
    166                ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
    167                ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
    168                IF(ZRIW.LT.GRCRIT) THEN
    169                  ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
    170                  ZB=1./GRCRIT+2./ZSQR
    171                  ZALPHA=0.5*(-ZB+SQRT(ZDEL))
    172                  ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK)
    173                  PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
    174                ELSE
    175                  PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
    176                ENDIF
    177             PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
    178           ENDIF
    179           ENDIF
    180   431 CONTINUE
    181  
    182   430 CONTINUE
    183136  440 CONTINUE
    184137 
    185 c     write(*,*) 'ptau'
    186 c     write(*,99) ((ji,ilevh,ptau(ji,ilevh),ji=1,NDLO2),
    187 c    .                  ilevh=1,klev+1)
    188  99   FORMAT(i3,i3,f15.5)
     138!     write(*,*) 'ZTAU'
     139!     write(*,99) ((ji,ilevh,ZTAU(ji,ilevh),ji=1,ndomainsz  ),
     140!    .                  ilevh=1,nlayer+1)
     141 99   FORMAT(i3,i3,f15.5)  ! not used
    189142
     143      !  Proganisation of the stress profile if breaking occurs at low level
     144      DO ji=1,kgwd
     145         jl=kdx(ji)
     146         Z_TAU(JL,IKENVH(JL))=ZTAU(JL,IKENVH(JL))
     147         Z_TAU(JL,IKCRITH(JL))=ZTAU(JL,IKCRITH(JL))
     148      ENDDO     
    190149
    191 C  REORGANISATION OF THE STRESS PROFILE
    192 C  IF BREAKING OCCURS AT LOW LEVEL:
    193 
    194       DO 530 ji=1,kgwd
    195       jl=kdx(ji)
    196       ZTAU(JL,KKENVH(JL))=PTAU(JL,KKENVH(JL))
    197       ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
    198  530  CONTINUE     
    199 
    200       DO 531 JK=1,KLEV
    201      
    202       DO 532 ji=1,kgwd
    203       jl=kdx(ji)
    204                
    205          IF(JK.GT.KKCRITH(JL).AND.JK.LT.KKENVH(JL))THEN
    206 
    207           ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKENVH(JL))
    208           ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KKENVH(JL))
    209           PTAU(JL,JK)=ZTAU(JL,KKENVH(JL)) +
    210      .                (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KKENVH(JL)) )*
    211      .                ZDELP/ZDELPT
    212      
    213         ENDIF
    214            
    215  532  CONTINUE   
    216  
    217  531  CONTINUE       
     150      DO JK=1,nlayer     
     151        DO ji=1,kgwd
     152           jl=kdx(ji)               
     153           IF(JK.GT.IKCRITH(JL).AND.JK.LT.IKENVH(JL))THEN
     154             ZDELP=pplev(JL,JK)-pplev(JL,IKENVH(JL))
     155             ZDELPT=pplev(JL,IKCRITH(JL))-pplev(JL,IKENVH(JL))
     156             ZTAU(JL,JK)=Z_TAU(JL,IKENVH(JL))+(Z_TAU(JL,IKCRITH(JL))-Z_TAU(JL,IKENVH(JL)))*ZDELP/ZDELPT   
     157          ENDIF           
     158          ENDDO
     159        end DO       
    218160
    219161      END SUBROUTINE GWPROFIL
    220162
    221       END MODULE gwprofil_mod
     163END MODULE gwprofil_mod

  • trunk/LMDZ.MARS/libf/phymars/gwstress_mod.F90

    r2641 r2642  
    1       MODULE gwstress_mod
     1MODULE gwstress_mod
    22     
    3       IMPLICIT NONE
     3IMPLICIT NONE
    44     
    5       CONTAINS
     5CONTAINS
    66
    7       SUBROUTINE GWSTRESS
    8      *         (  klon  , klev
    9      *         , KKCRIT, KSECT, KKHLIM, KTEST, KKCRITH, KCRIT, kkenvh
    10      *         , kknu
    11      *         , PRHO  , PSTAB, PVPH  , PVAR ,PVARor, psig 
    12      *         , PTFR  , PTAU 
    13      *         ,pgeom1 , pgamma, pd1  , pd2   ,pdmod ,pnu )
    14 C
    15 C**** *GWSTRESS*
    16 C
    17 C     PURPOSE.
    18 C     --------
    19 C
    20 C**   INTERFACE.
    21 C     ----------
    22 C     CALL *GWSTRESS*  FROM *GWDRAG*
    23 C
    24 C        EXPLICIT ARGUMENTS :
    25 C        --------------------
    26 C     ==== INPUTS ===
    27 C     ==== OUTPUTS ===
    28 C
    29 C        IMPLICIT ARGUMENTS :   NONE
    30 C        --------------------
    31 C
    32 C     METHOD.
    33 C     -------
    34 C
    35 C
    36 C     EXTERNALS.
    37 C     ----------
    38 C
    39 C
    40 C     REFERENCE.
    41 C     ----------
    42 C
    43 C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
    44 C
    45 C     AUTHOR.
    46 C     -------
    47 C
    48 C     MODIFICATIONS.
    49 C     --------------
    50 C     F. LOTT PUT THE NEW GWD ON IFS      22/11/93
    51 C
    52 C-----------------------------------------------------------------------
    53       use dimradmars_mod, only: ndlo2
     7      SUBROUTINE GWSTRESS(ngrid,nlayer,ktest,zrho, BV, pvar,psig,zgeom,zdmod,&
     8                ! Notice that this 3 variables are actually not used due to illness lines
     9                 ICRIT,IKNU,ZVPH,                                            &
     10                ! not used variables
     11 !                IKCRIT,ISECT,IKHLIM, IKCRITH,IKENVH,PVAR1,pgam,zd1,zd2,znu, &
     12                ! not defined not used variables
     13 !                 ZTFR
     14                !in(as 0.0)-output:
     15                ZTAU )
     16               
     17      !----------------------------------------------------------------------------------------------
     18      ! MODULE contains SUBROUTINE gwstress to compute low level stresses using subcritical, super
     19      ! critical forms.
     20      ! F. LOTT PUT THE NEW GWD ON IFS      22/11/93
     21      ! REFERENCE.
     22      ! SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
     23      ! Rewirten by J.Liu 03/03/2022
     24      !----------------------------------------------------------------------------------------------
     25
     26      use dimradmars_mod, only: ndomainsz   
    5427      implicit none
    55       integer klon,klev,kidia,kfdia
     28      include "yoegwd.h"
     29     
     30      ! 0. DECLARATIONS:
     31     
     32      ! 0.1   ARGUMENTS
     33      integer,intent(in):: ngrid    ! number of atmospheric columns
     34      integer,intent(in):: nlayer   ! number of atmospheric layers
     35      INTEGER,intent(in):: ktest(ndomainsz)   ! map of calling points
     36!      integer,intent(in):: IKCRIT(ndomainsz) ! not used
     37      integer,intent(in):: ICRIT(ndomainsz)   ! actually not used
     38 !     integer,intent(in):: IKCRITH(ndomainsz)! not used
     39 !     integer,intent(in):: ISECT(ndomainsz)  ! not used
     40 !     integer,intent(in):: IKHLIM(ndomainsz) ! not used
     41 !     integer,intent(in):: IKENVH(ndomainsz) ! The line use this variable has been commented
     42      integer,intent(in):: IKNU(ndomainsz)    ! actually not used
     43      REAL,INTENT(IN):: ZRHO(ndomainsz,nlayer+1) ! Density at 1/2 level
     44      REAL,INTENT(IN):: BV(ndomainsz,nlayer+1)   ! Brunt–Väisälä frequency at 1/2 level
     45      REAL,INTENT(IN):: ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H
     46      REAL,INTENT(IN):: zgeom(ndomainsz,nlayer)  ! Geopotetial height
     47      REAL,INTENT(IN):: pvar(ndomainsz)          ! Sub-grid scale standard deviation
     48!      REAL,INTENT(IN):: zd1(ndomainsz)       ! not used
     49!      REAL,INTENT(IN):: zd2(ndomainsz)       ! not used
     50!      REAL,INTENT(IN):: znu(ndomainsz)       ! not used
     51      REAL,INTENT(IN):: psig(ndomainsz)       ! SUB-GRID SCALE SLOPE
     52!      REAL,INTENT(IN):: pgam(ndomainsz)      ! not used
     53      REAL,INTENT(IN):: zdmod(ndomainsz)      ! Squre root of tao1 and tao2 without the constant, see equation 17 or 18
     54!      REAL,INTENT(IN):: ZTFR(ndomainsz)      ! not used. It is not even defined in this rountine
     55      REAL,INTENT(INOUT):: ZTAU(ndomainsz,nlayer+1) !GRAVITY WAVE STRESS.
     56     
     57      !0.2   LOCAL ARRAYS
     58      integer jl
     59      INTEGER kidia,kfdia
     60      real zvar    ! Sub-grid scale standard deviation at the calling points
     61      real zblock,zeff
     62      logical lo   ! actually not used bucause the if-endif condition that use this
     63                   ! variable has been commented
     64                   
     65!---------------------------------------------------------------------------------------------------
     66! 1. INITIALIZATION (not important initialization at all may be delete in the future)
     67!---------------------------------------------------------------------------------------------------
     68      kidia=1
     69      kfdia=ngrid
     70 100   CONTINUE ! continue tag without source, maybe need delete in future
     71!*         3.1     Gravity wave stress
     72 300   CONTINUE ! continue tag without source, maybe need delete in future
    5673
    57       include "yoegwd.h"
     74      DO JL=kidia,kfdia
     75        IF(KTEST(JL).EQ.1) THEN
     76        !Effective mountain height above the blocked flow     
     77!        IF(IKENVH(JL).EQ.nlayer)THEN
     78         ZBLOCK=0.0
     79!        ELSE
     80!         ZBLOCK=(zgeom(JL,IKENVH(JL))+zgeom(JL,IKENVH(JL)+1))/2./RG         
     81!        ENDIF   
     82        ZVAR=pvar(JL)
     83        ZEFF=AMAX1(0.,2.*ZVAR-ZBLOCK)
     84        ! Evaluate equation 17 to get the GW stress
     85        ZTAU(JL,nlayer+1)=zrho(JL,nlayer+1)*GKDRAG*psig(jl)*ZEFF**2    &
     86         /4./ZVAR*ZVPH(JL,nlayer+1)*zdmod(jl)*sqrt(BV(jl,nlayer+1))
    5887
    59 C-----------------------------------------------------------------------
    60 C
    61 C*       0.1   ARGUMENTS
    62 C              ---------
    63 C
    64       INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2),KSECT(NDLO2),
    65      *        KKHLIM(NDLO2),KTEST(NDLO2),KKENVH(NDLO2),KKNU(NDLO2)
    66 C
    67       REAL PRHO(NDLO2,klev+1),PSTAB(NDLO2,klev+1),PTAU(NDLO2,klev+1),
    68      *     PVPH(NDLO2,klev+1),PVAR(NDLO2,4),PTFR(NDLO2),
    69      *     pgeom1(NDLO2,klev),PVARor(NDLO2)
    70 C
    71       real pd1(NDLO2),pd2(NDLO2),pnu(NDLO2),psig(NDLO2),pgamma(NDLO2)
    72       real pdmod(NDLO2)
    73 C
    74 C-----------------------------------------------------------------------
    75 C
    76 C*       0.2   LOCAL ARRAYS
    77 C              ------------
    78       integer jl
    79       real zblock,zvar,zeff
    80       logical lo
    81 
    82 C
    83 C-----------------------------------------------------------------------
    84 C
    85 C*       0.3   FUNCTIONS
    86 C              ---------
    87 C     ------------------------------------------------------------------
    88 C
    89 C*         1.    INITIALIZATION
    90 C                --------------
    91 
    92 
    93       kidia=1
    94       kfdia=klon
    95 
    96 
    97 C
    98  100  CONTINUE
    99 C
    100 C*         3.1     GRAVITY WAVE STRESS.
    101 C
    102   300 CONTINUE
    103 C
    104 C
    105       DO 301 JL=kidia,kfdia
    106       IF(KTEST(JL).EQ.1) THEN
    107      
    108 C  EFFECTIVE MOUNTAIN HEIGHT ABOVE THE BLOCKED FLOW
    109  
    110 c        IF(KKENVH(JL).EQ.KLEV)THEN
    111          ZBLOCK=0.0
    112 c        ELSE
    113 c         ZBLOCK=(PGEOM1(JL,KKENVH(JL))+PGEOM1(JL,KKENVH(JL)+1))/2./RG         
    114 c        ENDIF
    115      
    116         ZVAR=PVAROR(JL)
    117         ZEFF=AMAX1(0.,2.*ZVAR-ZBLOCK)
    118 
    119         PTAU(JL,KLEV+1)=PRHO(JL,KLEV+1)*GKDRAG*psig(jl)*ZEFF**2
    120      *    /4./ZVAR*PVPH(JL,KLEV+1)*pdmod(jl)*sqrt(pstab(jl,klev+1))
    121 
    122 C  TOO SMALL VALUE OF STRESS OR  LOW LEVEL FLOW INCLUDE CRITICAL LEVEL
    123 C  OR LOW LEVEL FLOW:  GRAVITY WAVE STRESS NUL.
    124                
    125         LO=(PTAU(JL,KLEV+1).LT.GTSEC).OR.(KCRIT(JL).GE.KKNU(JL))
    126      *      .OR.(PVPH(JL,KLEV+1).LT.GVCRIT)
    127 c       IF(LO) PTAU(JL,KLEV+1)=0.0
    128      
    129       ELSE
    130      
    131           PTAU(JL,KLEV+1)=0.0
    132          
    133       ENDIF
    134      
    135   301 CONTINUE
    136 C
     88      !  Too small value of stress or low level flow include critical level
     89      !  or low level flow: gravity wave stress nul.                 
     90        LO=(ZTAU(JL,nlayer+1).LT.GTSEC).OR.(ICRIT(JL).GE.IKNU(JL)).OR. &
     91        (ZVPH(JL,nlayer+1).LT.GVCRIT)
     92!       IF(LO) ZTAU(JL,nlayer+1)=0.0     
     93        ELSE     
     94          ZTAU(JL,nlayer+1)=0.0         
     95        ENDIF     
     96      ENDDO
    13797
    13898      END SUBROUTINE GWSTRESS
    13999     
    140       END MODULE gwstress_mod
     100END MODULE gwstress_mod

  • trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F90

    r2641 r2642  
    1       MODULE orodrag_mod
    2      
    3       IMPLICIT NONE
    4      
    5       CONTAINS
    6      
    7       SUBROUTINE ORODRAG( klon,klev
    8      I                 , KGWD, KGWDIM, KDX, KTEST
    9      R                 , PTSPHY
    10      R                 , PAPHM1,PAPM1,PGEOM1,PTM1,PUM1
    11      R                 , PVM1, PVAROR, PSIG, PGAMMA, PTHETA
    12 C OUTPUTS
    13      R                 , PULOW,PVLOW
    14      R                 , PVOM,PVOL,PTE )
    15 C
    16 C
    17 C**** *ORODRAG* - DOES THE GRAVITY WAVE PARAMETRIZATION.
    18 C
    19 C     PURPOSE.
    20 C     --------
    21 C
    22 C          THIS ROUTINE COMPUTES THE PHYSICAL TENDENCIES OF THE
    23 C          PROGNOSTIC VARIABLES U,V  AND T DUE TO  VERTICAL TRANSPORTS BY
    24 C          SUBGRIDSCALE OROGRAPHICALLY EXCITED GRAVITY WAVES
    25 C
    26 C     EXPLICIT ARGUMENTS :
    27 C     --------------------
    28 C
    29 C     INPUT :
    30 C
    31 C     NLON               : NUMBER OF HORIZONTAL GRID POINTS
    32 C     NLEV               : NUMBER OF LEVELS
    33 C     KGWD               : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
    34 C     KGWDIM             : NUMBER OF POINTS AT WHICH THE SCHEME IS CALLED
    35 C     KDX(NLON)          : POINTS AT WHICH TO CALL THE SCHEME
    36 C     KTEST(NLON)        : MAP OF CALLING POINTS
    37 C     PTSPHY             : LENGTH OF TIME STEP
    38 C     PAPHM1(NLON,NLEV+1): PRESSURE AT MIDDLE LEVELS
    39 C     PAPM1(NLON,NLEV)   : PRESSURE ON MODEL LEVELS
    40 C     PGEOM1(NLON,NLEV)  : GEOPOTENTIAL HIEGHT OF MODEL LEVELS
    41 C     PTM1(NLON,NLEV)    : TEMPERATURE
    42 C     PUM1(NLON,NLEV)    : ZONAL WIND
    43 C     PVM1(NLON,NLEV)    : MERIDIONAL WIND
    44 C     PVAROR(NLON)       : SUB-GRID SCALE STANDARD DEVIATION
    45 C     PSIG(NLON)         : SUB-GRID SCALE SLOPE
    46 C     PGAMMA(NLON)       : SUB-GRID SCALE ANISOTROPY
    47 C     PTHETA(NLON)       : SUB-GRID SCALE PRINCIPAL AXES ANGLE
    48 C
    49 C     OUTPUT :
    50 C
    51 C     PULOW(NLON)        : LOW LEVEL ZONAL WIND
    52 C     PVLOW(NLON)        : LOW LEVEL MERIDIONAL WIND
    53 C     PVOM(NLON,NLEV)    : ZONAL WIND TENDENCY
    54 C     PVOL(NLON,NLEV)    : MERIDIONAL WIND TENDENCY
    55 C     PTE(NLON,NLEV)     : TEMPERATURE TENDENCY
    56 C
    57 C     IMPLICIT ARGUMENTS :
    58 C     --------------------
    59 C
    60 C     comcstfi.h
    61 C     yoegwd.h
    62 C
    63 C     METHOD.
    64 C     -------
    65 C
    66 C     EXTERNALS.
    67 C     ----------
    68 C
    69 C     REFERENCE.
    70 C     ----------
    71 C
    72 C     AUTHOR.
    73 C     -------
    74 C     M.MILLER + B.RITTER   E.C.M.W.F.     15/06/86.
    75 C
    76 C     F.LOTT + M. MILLER    E.C.M.W.F.     22/11/94
    77 C-----------------------------------------------------------------------
    78       use dimradmars_mod, only: ndlo2
     1MODULE orodrag_mod
     2     
     3IMPLICIT NONE
     4     
     5CONTAINS
     6     
     7      SUBROUTINE ORODRAG( ngrid,nlayer, kgwd, kgwdim, kdx, ktest, ptimestep, &
     8                  pplev,pplay,zgeom,pt,pu, pv, pvar, psig, pgam, pthe,       &
     9                  !OUTPUTS
     10                  PULOW,PVLOW, pdudt,pdvdt,pdtdt)
     11                 
     12     !--------------------------------------------------------------------------------
     13     ! The main routine ORODRAG that does the orographic gravity wave parameterization.
     14     ! This routine computes the physical tendencies of the prognostic variables U,V, and
     15     ! T due to vertical transport by subgridscale orographically excited gravity waves
     16     ! M.MILLER + B.RITTER   E.C.M.W.F.     15/06/86.
     17     ! F.LOTT + M. MILLER    E.C.M.W.F.     22/11/94
     18     !--
     19     ! The purpose of this routine is:
     20     ! 1) call OROSETUP to get the parameters such as the top altitude(level) of the blocked
     21     !    flow and other critical levels.
     22     ! 2) call GWSTRESS and GWPROFILE to get the gw stress.
     23     ! 3) calculate the zonal/meridional wind tendency based on the GW stress and mountain drag
     24     !    on the flow.
     25     ! 4) calculate the temperature tendency based on the differtial of square of the wind
     26     !    between t and t+dt.
     27     ! Rewrited by J.LIU     LMDZ.COMMON/phymars/libf  29/01/2022
     28     !--
     29     ! REFERENCE.
     30     ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
     31     ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization:
     32     !    Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society,
     33     !    123(537), 101-127.
     34     !-----------------------------------------------------------------------
     35     
     36      use dimradmars_mod, only: ndomainsz
    7937      USE gwstress_mod, ONLY: gwstress
    8038      USE gwprofil_mod, ONLY: gwprofil
    8139      USE comcstfi_h, ONLY: g, cpp
     40     
    8241      implicit none
    83 C
    84 C
    85       integer klon,klev,kidia
    86       parameter(kidia=1)
    87       integer, save :: kfdia ! =NDLO2
    88      
    89 !$OMP THREADPRIVATE(kfdia)
    90 
    9142      include "yoegwd.h"
    92 C-----------------------------------------------------------------------
    93 C
    94 C*       0.1   ARGUMENTS
    95 C              ---------
    96 C
    97 C
    98       REAL  PTE(NDLO2,klev),
    99      *      PVOL(NDLO2,klev),
    100      *      PVOM(NDLO2,klev),
    101      *      PULOW(NDLO2),
    102      *      PVLOW(NDLO2)
    103       REAL  PUM1(NDLO2,klev),
    104      *      PVM1(NDLO2,klev),
    105      *      PTM1(NDLO2,klev),
    106      *      PVAROR(NDLO2),PSIG(NDLO2),PGAMMA(NDLO2),PTHETA(NDLO2),
    107      *      PGEOM1(NDLO2,klev),
    108      *      PAPM1(NDLO2,klev),
    109      *      PAPHM1(NDLO2,klev+1)
    110 C
    111       integer kgwd,kgwdim
    112       real ptsphy
    113       INTEGER  KDX(NDLO2),KTEST(NDLO2)
    114 C-----------------------------------------------------------------------
    115 C
    116 C*       0.2   LOCAL ARRAYS
    117 C              ------------
    118       INTEGER  ISECT(NDLO2),
    119      *         ICRIT(NDLO2),
    120      *         IKCRITH(NDLO2),
    121      *         IKenvh(NDLO2),
    122      *         IKNU(NDLO2),
    123      *         IKNU2(NDLO2),
    124      *         IKCRIT(NDLO2),
    125      *         IKHLIM(NDLO2)
    126       integer ji,jk,jl,klevm1,ilevp1
    127 C      real gkwake
    128       real ztmst,pvar(NDLO2,4),ztauf(NDLO2,klev+1)
    129       real zrtmst,zdelp,zb,zc,zbet
    130       real zconb,zabsv,zzd1,ratio,zust,zvst,zdis,ztemp
    131 C
    132       REAL   ZTAU(NDLO2,klev+1),
    133      *       ZSTAB(NDLO2,klev+1),
    134      *       ZVPH(NDLO2,klev+1),
    135      *       ZRHO(NDLO2,klev+1),
    136      *       ZRI(NDLO2,klev+1),
    137      *       ZpsI(NDLO2,klev+1),
    138      *       Zzdep(NDLO2,klev)
    139       REAL   ZDUDT(NDLO2),
    140      *       ZDVDT(NDLO2),
    141      *       ZDTDT(NDLO2),
    142      *       ZDEDT(NDLO2),
    143      *       ZVIDIS(NDLO2),
    144      *       ZTFR(NDLO2),
    145      *       Znu(NDLO2),
    146      *       Zd1(NDLO2),
    147      *       Zd2(NDLO2),
    148      *       Zdmod(NDLO2)
    149 C
    150 C------------------------------------------------------------------
    151 C
    152 C*         1.    INITIALIZATION
    153 C                --------------
    154 C
    155  100  CONTINUE
    156 C
    157 C     ------------------------------------------------------------------
    158 C
    159 C*         1.1   COMPUTATIONAL CONSTANTS
    160 C                -----------------------
    161 C
    162  110  CONTINUE
    163 C
    164       kfdia=NDLO2
    165      
    166 c     ZTMST=TWODT
    167 c     IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT
    168       KLEVM1=KLEV-1
    169       ZTMST=PTSPHY
    170       ZRTMST=1./ZTMST
    171 C     ------------------------------------------------------------------
    172 C
    173  120  CONTINUE
    174 C
    175 C     ------------------------------------------------------------------
    176 C
    177 C*         1.3   CHECK WHETHER ROW CONTAINS POINT FOR PRINTING
    178 C                ---------------------------------------------
    179 C
    180  130  CONTINUE
    181 C
    182 C     ------------------------------------------------------------------
    183 C
    184 C*         2.     PRECOMPUTE BASIC STATE VARIABLES.
    185 C*                ---------- ----- ----- ----------
    186 C*                DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
    187 C*                LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
    188 C*                THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
    189 C
    190   200 CONTINUE
    191 C
    192 C
    193 C
    194       CALL OROSETUP
    195      *     ( klon, klev , KTEST
    196      *     , IKCRIT, IKCRITH, ICRIT, ISECT, IKHLIM, ikenvh,IKNU,iknu2
    197      *     , PAPHM1, PAPM1 , PUM1   , PVM1 , PTM1 , PGEOM1, pvaror
    198      *     , ZRHO  , ZRI   , ZSTAB  , ZTAU , ZVPH , zpsi, zzdep
    199      *     , PULOW, PVLOW
    200      *     , ptheta,pgamma,znu  ,zd1,  zd2,  zdmod                )
    201 C
    202 C
    203 C
    204 C***********************************************************
    205 C
    206 C
    207 C*         3.      COMPUTE LOW LEVEL STRESSES USING SUBCRITICAL AND
    208 C*                 SUPERCRITICAL FORMS.COMPUTES ANISOTROPY COEFFICIENT
    209 C*                 AS MEASURE OF OROGRAPHIC TWODIMENSIONALITY.
    210 C
    211   300 CONTINUE
    212 C
    213       CALL GWSTRESS
    214      *    ( klon  , klev
    215      *    , IKCRIT, ISECT, IKHLIM, KTEST, IKCRITH, ICRIT, ikenvh, IKNU
    216      *    , ZRHO  , ZSTAB, ZVPH  , PVAR , pvaror,  psig
    217      *    , ZTFR   , ZTAU
    218      *    , pgeom1,pgamma,zd1,zd2,zdmod,znu)
    219 C
    220 C*         4.      COMPUTE STRESS PROFILE.
    221 C*                 ------- ------ --------
    222 C
    223   400 CONTINUE
    224 C
    225 C
    226       CALL GWPROFIL
    227      *       (  klon , klev
    228      *       , kgwd   , kdx  , KTEST
    229      *       , IKCRIT, IKCRITH, ICRIT  , ikenvh, IKNU
    230      *       ,iknu2 , pAPHM1, ZRHO   , ZSTAB , ZTFR   , ZVPH
    231      *       , ZRI   , ZTAU   , ztauf
    232      *       , zdmod , znu    , psig  , pgamma , pvaror   )
    233 C
    234 C
    235 C*         5.      COMPUTE TENDENCIES.
    236 C*                 -------------------
    237 C
    238   500 CONTINUE
    239 C
    240 C  EXPLICIT SOLUTION AT ALL LEVELS FOR THE GRAVITY WAVE
    241 C  IMPLICIT SOLUTION FOR THE BLOCKED LEVELS
    242 
    243       DO 510 JL=KIDIA,KFDIA
    244       ZVIDIS(JL)=0.0
    245       ZDUDT(JL)=0.0
    246       ZDVDT(JL)=0.0
    247       ZDTDT(JL)=0.0
    248   510 CONTINUE
    249 C
    250       ILEVP1=KLEV+1
    251 C
    252 C
    253       DO 524 JK=1,klev
    254 C
    255 CDIR$ IVDEP
    256 C
    257 C      GKWAKE=0.5
    258 C
    259 C     NOW SET IN SUGWD.F
    260 C
    261       DO 523 JL=1,KGWD
    262       JI=KDX(JL)
    263       ZDELP=pAPHM1(Ji,JK+1)-pAPHM1(Ji,JK)
    264       ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP)
    265       ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
    266       ZDVDT(JI)=(pvLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
    267       if(jk.ge.ikenvh(ji)) then
    268          zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2
    269          zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2
    270          zconb=2.*ztmst*GKWAKE*psig(ji)/(4.*pvaror(ji))
    271          zabsv=sqrt(PUM1(JI,JK)**2+PVM1(JI,JK)**2)/2.
    272          zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
    273          ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/
    274      *   (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
    275          zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
    276          zdudt(ji)=-pum1(ji,jk)/ztmst
    277          zdvdt(ji)=-pvm1(ji,jk)/ztmst
    278          zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
    279          zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
    280       end if
    281       PVOM(JI,JK)=ZDUDT(JI)
    282       PVOL(JI,JK)=ZDVDT(JI)
    283       ZUST=PUM1(JI,JK)+ZTMST*ZDUDT(JI)
    284       ZVST=PVM1(JI,JK)+ZTMST*ZDVDT(JI)
    285       ZDIS=0.5*(PUM1(JI,JK)**2+PVM1(JI,JK)**2-ZUST**2-ZVST**2)
    286       ZDEDT(JI)=ZDIS/ZTMST
    287       ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP
    288       ZDTDT(JI)=ZDEDT(JI)/cpp
    289       PTE(JI,JK)=ZDTDT(JI)
    290 
    291   523 CONTINUE
    292 
    293   524 CONTINUE
    294 C
    295 C
     43     
     44      ! 0. DECLARATIONS:
     45     
     46      ! 0.1   INPUTS
     47      INTEGER, intent(in):: ngrid          ! Number of atmospheric columns
     48      INTEGER, intent(in):: nlayer         ! Number of atmospheric layers
     49      INTEGER, intent(in):: kgwd           ! Number of points at which the scheme is called
     50      INTEGER, intent(in):: kgwdim         ! kgwdim=max(1,kgwd)
     51      INTEGER, intent(in):: kdx(ndomainsz)     ! Points at which to call the scheme.
     52      INTEGER, intent(in):: ktest(ndomainsz)   ! Map of calling points   
     53
     54      REAL, intent(in):: pu(ndomainsz,nlayer)  ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu)
     55      REAL, intent(in):: pv(ndomainsz,nlayer)  ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv)
     56      REAL, intent(in):: pt(ndomainsz,nlayer)  ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt)
     57      REAL, intent(in):: pvar(ndomainsz)       ! Sub-grid scale standard deviation
     58      REAL, intent(in):: psig(ndomainsz)       ! Sub-grid scale slope
     59      REAL, intent(in):: pgam(ndomainsz)       ! Sub-grid scale anisotropy
     60      REAL, intent(in):: pthe(ndomainsz)       ! Sub-grid scale principal axes angle
     61      REAL, intent(in):: zgeom(ndomainsz,nlayer)     ! Geopotential height of full levels
     62      REAL, intent(in):: pplay(ndomainsz,nlayer)     ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay)
     63      REAL, intent(in):: pplev(ndomainsz,nlayer+1)   ! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev) 
     64      REAL, intent(in):: ptimestep                   ! Time step of the Physics(s)
     65     
     66      ! 0.2   OUTPUTS
     67      REAL, intent(out):: pdtdt(ndomainsz,nlayer)   ! Tendency on temperature (K/s) due to orographic gravity waves
     68      REAL, intent(out):: pdvdt(ndomainsz,nlayer)   ! Tendency on meridional wind (m/s/s) due to orographic gravity waves
     69      REAL, intent(out):: pdudt(ndomainsz,nlayer)   ! Tendency on zonal wind (m/s/s) due to orographic gravity waves
     70      REAL, intent(out):: PULOW(ndomainsz)          ! Low level zonal wind
     71      REAL, intent(out):: PVLOW(ndomainsz)          ! Low level meridional wind
     72
     73      !0.3   LOCAL ARRAYS             
     74!      INTEGER ISECT(ndomainsz)      ! not used ?
     75      INTEGER ICRIT(ndomainsz)       ! Critical layer where orographic GW breaks
     76      INTEGER IKCRITH(ndomainsz)     ! Dynamical mixing height for the breaking of gravity waves
     77      INTEGER IKenvh(ndomainsz)      ! Top of the  blocked layer
     78      INTEGER IKNU(ndomainsz)        ! 4*pvar layer
     79      INTEGER IKNU2(ndomainsz)       ! 3*pvar layer
     80      INTEGER IKCRIT(ndomainsz)      ! Top of low level flow height
     81!      INTEGER IKHLIM(ndomainsz)     ! not used?
     82     
     83      integer ji,jk,jl
     84      integer ilevp1                    !=nlayer+1 just used once. Maybe directly use nlayer+1 to replace in the future
     85      
     86!      real ztauf(ndomainsz,nlayer+1)   ! not used
     87      real zdelp                        ! =pplev(nlayer+1)-pplev(nlayer) dp of two 1/2 levels
     88      real zb                           ! B(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2)
     89      real zc                           ! C(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2)
     90      real zbet                         ! Equation (16) blocked flow drag
     91      real zconb                        ! Middle part of the equation (16)
     92      real zabsv                        ! Tail of equation (16),U*|U|/2+V*|V|/2
     93      real zzd1                         ! Second tail part of the equation (16)
     94      real ratio                        ! Aspect ratio of the obstacle (mountain),equation (14)
     95      real zust                         ! U(t+dt)
     96      real zvst                         ! V(t+dt)
     97      real zdis                         ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|]
     98      real ztemp                        ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress?
     99      REAL ZTAU(ndomainsz,nlayer+1)   ! Output of subrountine GWPROFILE: Gravtiy wave stress
     100      REAL BV(ndomainsz,nlayer+1)     ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL)
     101      REAL ZVPH(ndomainsz,nlayer+1)   ! Low level wind speed U_H in Ref.2
     102      REAL ZRHO(ndomainsz,nlayer+1)   ! density calculated (output) by the OROSETUP rountines(but not used by ORODRAG)
     103                                      ! used as input for GWSTRESS and GWPROFIL
     104      REAL PRI(ndomainsz,nlayer+1)    ! Mean flow richardson number at 1/2 level
     105      REAL ZpsI(ndomainsz,nlayer+1)   ! The angle between the incident flow direction and the normal ridge direction of the mountain
     106      REAL Zzdep(ndomainsz,nlayer)    ! dp by full level
     107      REAL ZDUDT(ndomainsz)           ! du/dt due to oro-gw
     108      REAL ZDVDT(ndomainsz)           ! dv/dt due to oro-gw
     109      REAL ZDTDT(ndomainsz)           ! dT/dt due to oro-gw
     110      REAL ZDEDT(ndomainsz)           ! zdis/ptimestemp, zdedt/Cp=dT/dt
     111      REAL ZVIDIS(ndomainsz)          ! in an invalid line, not used
     112      REAL Znu(ndomainsz)             ! A critical value see equation 9(Since it only compute inside OROSETUP,Maybe needs delete in future)
     113      REAL Zd1(ndomainsz)             ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18 in Ref.2
     114      REAL Zd2(ndomainsz)             ! (B-C)sin(psi)cos(psi)   see equation 17 or 18 in Ref.2
     115      REAL Zdmod(ndomainsz)           ! sqrt(zd1^2+zd2^2)
     116
     117!------------------------------------------------------------------------------------
     118! 1.INITIALIZATION (NOT USED )
     119!------------------------------------------------------------------------------------
     120! 100  CONTINUE   ! continue tag without source, maybe need delete in future
     121!*         1.1   Computional constants
     122! 110  CONTINUE ! continue tag without source, maybe need delete in future
     123!      kfdia=ndomainsz     
     124!     ZTMST=TWODT
     125!     IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT
     126!      nlayerM1=nlayer-1
     127!      ZTMST=ptimestep
     128!      ZRTMST=1./ptimestep
     129! 120  CONTINUE ! continue tag without source, maybe need delete in future
     130!*         1.3   Check whether row contains point for printing
     131! 130  CONTINUE ! continue tag without source, maybe need delete in future
     132 
     133!------------------------------------------------------------------------------------
     134! 2. Precompute basic state variables .
     135!------------------------------------------------------------------------------------     
     136! 200  CONTINUE ! continue tag without source, maybe need delete in future
     137      ! Define low level wind,project winds in plane of low level wind, determine sector
     138      ! in which to take the variance and set indicator for critical levels       
     139      CALL OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom,       &
     140                           pvar,pthe, pgam,                                       &
     141                           !output in capital
     142                           IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2,             &
     143                           !ISECT, IKHLIM, not used
     144                           ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD,    &
     145                           PULOW, PVLOW)
     146
     147!------------------------------------------------------------------------------------
     148! 3. Computes low level stresses using subcritical and super critical forms.
     149!    Computes anisotropy coefficient as measure of orographic two-dimensionality
     150!------------------------------------------------------------------------------------           
     151!  300 CONTINUE ! continue tag without source, maybe need delete in future
     152      ! Computes low level stresses
     153      CALL GWSTRESS(ngrid,nlayer,ktest,zrho, BV, pvar,psig,zgeom,zdmod,&
     154                ! Notice that this 3 variables are actually not used due to illness lines
     155                 ICRIT,IKNU,ZVPH,                                      &
     156                ! not used variables
     157 !                IKCRIT,ISECT,IKHLIM, IKCRITH,IKENVH,PVAR1,pgam,zd1,zd2,znu, &
     158                ! not defined not used variables
     159 !                 ZTFR
     160                !in(as 0.0)-output:
     161                ZTAU )
     162               
     163!------------------------------------------------------------------------------------
     164! 4. Compute stress profile.
     165!------------------------------------------------------------------------------------
     166!  400 CONTINUE ! continue tag without source, maybe need delete in future
     167      ! Compute stress profile.
     168      CALL GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, &
     169                 IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev,  &
     170                 ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar,    &
     171                 !not used variables
     172                 !IKCRIT,ZTAUf,ZTFR,znu,pgam,
     173                 ! in-output
     174                 ZTAU      )
     175
     176!------------------------------------------------------------------------------------
     177! 5. Compute tendencies.
     178!------------------------------------------------------------------------------------
     179!  500 CONTINUE ! continue tag without source, maybe need delete in future
     180!  EXplicit solution at all levels for gravity wave
     181!  Implicit solution for the blocked levels
     182!      DO JL=KIDIA,KFDIA
     183      ! Initialization the tendencies
     184      DO JL=1,ndomainsz
     185            ZVIDIS(JL)=0.0  !An invalid lines contain this variable
     186            ZDUDT(JL)=0.0
     187            ZDVDT(JL)=0.0
     188            ZDTDT(JL)=0.0
     189      ENDDO
     190     
     191      ! all the calculation of equation (16-18) are in the flowing loop
     192      ! 1) if the wind flow over the mountain [Either the mountain is too low or the
     193      !    flow higher than the blocked level], then calculate the wind tendencies by oro-GW stress directly
     194      ! 2) if the flow level lower than the blocked level, then calculate the wind tendencies by obstacle drag
     195      !   
     196      ILEVP1=nlayer+1
     197      DO  JK=1,nlayer
     198        DO  JL=1,kgwd
     199             ! Pick up the position of each calling points
     200             JI=kdx(JL)
     201             ! dp
     202             ZDELP=pplev(Ji,JK+1)-pplev(Ji,JK)
     203             ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress?
     204             ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP)
     205             ! du/dt=(U_low*tao1-V_low*tao2)*T/(tao1^2+tao2^2)*0.5
     206             ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
     207             ! dv/dt=(V_low*tao1+U_low*tao2)*T/(tao1^2+tao2^2)*0.5
     208             ZDVDT(JI)=(PVLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
     209             if(jk.ge.ikenvh(ji)) then  ! IF the level lower than(notice here the level is inversed thus it is lower)
     210                                        ! the top of the blocked layer (commonlly is the highest top of a mountain.)
     211                                        ! Then use equation (16) to calculate the mountain drag on the incident flow
     212                ! Coefficient B and C Ref.2
     213                zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
     214                zc=0.48*pgam(ji)+0.3*pgam(ji)**2
     215                ! Middle part of the equation (16)
     216                zconb=2.*ptimestep*GKWAKE*psig(ji)/(4.*pvar(ji))
     217                ! Tail part of the equation (16),U*|U|/2+V*|V|/2
     218                zabsv=sqrt(pu(JI,JK)**2+pv(JI,JK)**2)/2.
     219                ! Second tail part of the equation (16)
     220                zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
     221                ! Aspect ratio of the obstacle(topography), equation(14)
     222                ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/(pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
     223                ! Equation (16)
     224                zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
     225                ! du/dt and dv/dt due to the mountain drag
     226                zdudt(ji)=-pu(ji,jk)/ptimestep
     227                zdvdt(ji)=-pv(ji,jk)/ptimestep
     228                zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
     229                zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
     230             end if
     231             ! wind tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked))
     232             pdudt(JI,JK)=ZDUDT(JI)
     233             pdvdt(JI,JK)=ZDVDT(JI)
     234             ! winds at t+dt (oro-gw induced increaments are added)
     235             ZUST=pu(JI,JK)+ptimestep*ZDUDT(JI)
     236             ZVST=pv(JI,JK)+ptimestep*ZDVDT(JI)
     237             ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|]
     238             ZDIS=0.5*(pu(JI,JK)**2+pv(JI,JK)**2-ZUST**2-ZVST**2)
     239             ZDEDT(JI)=ZDIS/ptimestep
     240             ! This line no use maybe delete in the future
     241             ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP
     242             ! Temperature tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked))
     243             ZDTDT(JI)=ZDEDT(JI)/cpp
     244             pdtdt(JI,JK)=ZDTDT(JI)
     245        ENDDO !JL=1,kgwd
     246      end DO !JK=1,nlayer
    296247
    297248      END SUBROUTINE ORODRAG
    298249     
    299       END MODULE orodrag_mod
     250END MODULE orodrag_mod

  • trunk/LMDZ.MARS/libf/phymars/orosetup.F90

    r2641 r2642  
    1       SUBROUTINE OROSETUP
    2      *         ( klon   , klev  , KTEST
    3      *         , KKCRIT, KKCRITH, KCRIT, KSECT , KKHLIM
    4      *         , kkenvh, kknu  , kknu2
    5      *         , PAPHM1, PAPM1 , PUM1   , PVM1 , PTM1  , PGEOM1, pvaror
    6      *         , PRHO  , PRI   , PSTAB  , PTAU , PVPH  ,ppsi, pzdep
    7      *         , PULOW , PVLOW 
    8      *         , Ptheta, pgamma,  pnu  ,  pd1  ,  pd2  ,pdmod  )
    9 C
    10 C**** *GWSETUP*
    11 C
    12 C     PURPOSE.
    13 C     --------
    14 C
    15 C**   INTERFACE.
    16 C     ----------
    17 C          FROM *ORODRAG*
    18 C
    19 C        EXPLICIT ARGUMENTS :
    20 C        --------------------
    21 C     ==== INPUTS ===
    22 C     ==== OUTPUTS ===
    23 C
    24 C        IMPLICIT ARGUMENTS :   NONE
    25 C        --------------------
    26 C
    27 C     METHOD.
    28 C     -------
    29 C
    30 C
    31 C     EXTERNALS.
    32 C     ----------
    33 C
    34 C
    35 C     REFERENCE.
    36 C     ----------
    37 C
    38 C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
    39 C
    40 C     AUTHOR.
    41 C     -------
    42 C
    43 C     MODIFICATIONS.
    44 C     --------------
    45 C     F.LOTT  FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993
    46 C
    47 C-----------------------------------------------------------------------
    48       use dimradmars_mod, only: ndlo2
     1SUBROUTINE OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom, &
     2            pvar,pthe, pgam,                                       &
     3            !output in capital
     4            IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2,             &
     5            !ISECT, IKHLIM, not used
     6            ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD,    &
     7            PULOW, PVLOW)
     8      !---------------------------------------------------------------------------------------------       
     9      !!**** *GWSETUP*! Computes low level stresses using subcritical and super critical forms.
     10      ! As well as, computes anisotropy coefficient as measure of orographic two-dimensionality
     11      ! F.LOTT  FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993!
     12      !--
     13      ! REFERENCE.
     14      ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
     15      ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization:
     16      !    Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society,
     17      !    123(537), 101-127.
     18      !--
     19      ! MODIFICATIONS.
     20      ! 1.Rewiten by J.liu 03/03/2022
     21      !-----------------------------------------------------------------------
     22      use dimradmars_mod, only: ndomainsz
    4923      USE comcstfi_h
    5024      implicit none
    51 C
    52 
    53       integer klon,klev,kidia,kfdia
    54 
    55 #include "yoegwd.h"
    56 
    57 C-----------------------------------------------------------------------
    58 C
    59 C*       0.1   ARGUMENTS
    60 C              ---------
    61 C
    62       INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2),KSECT(NDLO2),
    63      *        KKHLIM(NDLO2),KTEST(NDLO2),kkenvh(NDLO2)
    64 
    65 C
    66       REAL PAPHM1(NDLO2,KLEV+1),PAPM1(NDLO2,KLEV),PUM1(NDLO2,KLEV),
    67      *     PVM1(NDLO2,KLEV),PTM1(NDLO2,KLEV),PGEOM1(NDLO2,KLEV),
    68      *     PRHO(NDLO2,KLEV+1),PRI(NDLO2,KLEV+1),PSTAB(NDLO2,KLEV+1),
    69      *     PTAU(NDLO2,KLEV+1),PVPH(NDLO2,KLEV+1),ppsi(NDLO2,klev+1),
    70      *     pzdep(NDLO2,klev)
    71        REAL PULOW(NDLO2),PVLOW(NDLO2),ptheta(NDLO2),pgamma(NDLO2),
    72      *     pnu(NDLO2),
    73      *     pd1(NDLO2),pd2(NDLO2),pdmod(NDLO2)
    74       real pvaror(NDLO2)
    75 C
    76 C-----------------------------------------------------------------------
    77 C
    78 C*       0.2   LOCAL ARRAYS
    79 C              ------------
    80 C
    81 C
    82       LOGICAL LL1(NDLO2,klev+1)
    83       integer kknu(NDLO2),kknu2(NDLO2),kknub(NDLO2),kknul(NDLO2),
    84      *        kentp(NDLO2),ncount(NDLO2) 
    85 C
    86       REAL ZHCRIT(NDLO2,klev),ZNCRIT(NDLO2,klev),
    87      *     ZVPF(NDLO2,klev), ZDP(NDLO2,klev)
    88       REAL ZNORM(NDLO2),zpsi(NDLO2),zb(NDLO2),zc(NDLO2),
    89      *      zulow(NDLO2),zvlow(NDLO2),znup(NDLO2),znum(NDLO2)
    90 C
    91 c   declarations pour "implicit none"
    92       integer jk,jl,ilevm1,ilevm2,ilevh
    93       real zu,zphi,zcons1,zcons2,zcons3,zwind,zdwind,zhgeo
    94       real zvt1,zvt2,zst,zvar,zdelp,zstabm,zstabp,zrhom,zrhop
    95       real alpha,zggeenv,zggeom1,zgvar
     25      include "yoegwd.h"
     26     
     27      ! 0. DECLARATIONS:
     28     
     29      ! 0.1 inputs:
     30      integer,intent(in):: ngrid    ! number of atmospheric columns
     31      integer,intent(in):: nlayer   ! number of atmospheric layers
     32      INTEGER,intent(in):: ktest(ndomainsz)  ! map of calling points
     33
     34      real, intent(in) :: pplev(ndomainsz,nlayer+1)! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev)
     35      real, intent(in) :: pplay(ndomainsz,nlayer)  ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay)
     36      real, intent(in) :: pu(ndomainsz,nlayer)     ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu)
     37      real, intent(in) :: pv(ndomainsz,nlayer)     ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv)
     38      real, intent(in) :: pt(ndomainsz,nlayer)     ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt)
     39      real, intent(in) :: zgeom(ndomainsz,nlayer)  ! Geopotetial height
     40      real, intent(in) :: pvar(ndomainsz)          ! Sub-grid scale standard deviation
     41      real, intent(in) :: pthe(ndomainsz)          ! Sub-grid scale principal axes angle
     42      real, intent(inout) :: pgam(ndomainsz)       ! Sub-grid scale anisotropy
     43     
     44      ! 0.2 outputs:
     45      INTEGER,intent(out):: IKCRIT(ndomainsz)       ! top of low level flow height
     46      INTEGER,intent(out):: IKCRITH(ndomainsz)      ! dynamical mixing height for the breaking of gravity waves
     47      INTEGER,intent(out):: ICRIT(ndomainsz)        ! Critical layer where orographic GW breaks
     48!      INTEGER,intent(out):: ISECT(ndomainsz)       ! not used
     49!      INTEGER,intent(out):: IKHLIM(ndomainsz)      ! not used
     50      INTEGER,intent(out):: IKENVH(ndomainsz)       ! Top of the  blocked layer
     51      INTEGER,intent(out):: IKNU(ndomainsz)         ! 4*pvar layer
     52      INTEGER,intent(out):: IKNU2(ndomainsz)        ! 3*pvar layer
     53
     54      REAL, intent(out):: ZRHO(ndomainsz,nlayer+1)  ! Density at 1/2 level
     55      REAL, intent(out):: PRI(ndomainsz,nlayer+1)   ! Mean flow richardson number at 1/2 level
     56      REAL, intent(out):: BV(ndomainsz,nlayer+1)    ! Brunt–Väisälä frequency at 1/2 level
     57      REAL, intent(out):: ZTAU(ndomainsz,nlayer+1)  ! Gravity wave stress. Set to 0.0 here and will calculate in GWSTRESS later
     58      REAL, intent(out):: ZVPH(ndomainsz,nlayer+1)  ! Low level wind speed U_H
     59      REAL, intent(out):: ZPSI(ndomainsz,nlayer+1)  ! The angle between the incident flow direction and the normal ridge direction pthe
     60      REAL, intent(out):: ZZDEP(ndomainsz,nlayer)   ! dp by full level
     61     
     62      REAL, intent(out):: PULOW(ndomainsz)          ! Low level zonal wind
     63      REAL, intent(out):: PVLOW(ndomainsz)          ! Low level meridional wind
     64      REAL, intent(out):: ZNU(ndomainsz)            ! A critical value see equation 9
     65      REAL, intent(out):: ZD1(ndomainsz)            ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18
     66      REAL, intent(out):: ZD2(ndomainsz)            ! (B-C)sin(psi)cos(psi)   see equation 17 or 18
     67      REAL, intent(out):: ZDMOD(ndomainsz)          ! sqrt(zd1^2+zd2^2)
     68
     69      !0.3 Local arrays   
     70      integer IKNUb(ndomainsz)   ! 2*pvar layer
     71      integer IKNUl(ndomainsz)   ! 1*pvar layer
     72      integer kentp(ndomainsz)   ! initialized value but never used 
     73      integer ncount(ndomainsz)  ! initialized value but never used
     74
     75      REAL ZHCRIT(ndomainsz,nlayer)  ! tag for 1*pvar, 2*pvar,3*pvar and 4*pvar, pvar is mu means SD
     76!      REAL ZNCRIT(ndomainsz,nlayer) ! not used
     77      REAL ZVPF(ndomainsz,nlayer)    ! Flow in plane of low level stress. Seems a unit vector
     78      REAL ZDP(ndomainsz,nlayer)     ! dp differitial of pressure
     79     
     80      REAL ZNORM(ndomainsz)     ! The norm ridge of a moutain?
     81      REAL zb(ndomainsz)        ! Parameter B in eqution 17
     82      REAL zc(ndomainsz)        ! Parameter C in eqution 17
     83      REAL zulow(ndomainsz)     ! initialized value but never used     
     84      REAL zvlow(ndomainsz)     ! initialized value but never used 
     85      REAL znup(ndomainsz)      ! znu in top of 1/2 level
     86      REAL znum(ndomainsz)      ! znu in bottom of 1/2 level
     87
     88      integer jk,jl
     89      integer ilevm1 !=nlayer-1
     90      integer ilevm2 !=nlayer-2
     91      integer ilevh  !=nalyer/3
     92      INTEGER kidia  !=1
     93      INTEGER kfdia  !=ngrid
     94      real zcons1    !=1/r
     95      real zcons2    !=g^2/cpp
     96      real zcons3    !=1.5*pi
     97      real zphi      ! direction of the incident flow
     98      real zhgeo     ! Height calculated by geopotential/g
     99      real zu        ! Low level zonal wind (to denfine the dirction of background wind)
     100      real zwind,zdwind
     101      real zvt1,zvt2,zst,zvar
     102      real zdelp     !dp differitial of pressure
     103      ! variables for bv and density rho
     104      real zstabm,zstabp,zrhom,zrhop
     105!      real alpha    !=3. but never used
     106      real zggeenv,zggeom1,zgvar
    96107      logical lo
    97 
    98 C     ------------------------------------------------------------------
    99 C
    100 C*         1.    INITIALIZATION
    101 C                --------------
    102 C
    103 c     print *,' entree gwsetup'
    104  100  CONTINUE
    105 C
    106 C     ------------------------------------------------------------------
    107 C
    108 C*         1.1   COMPUTATIONAL CONSTANTS
    109 C                -----------------------
    110 C
    111 
     108      LOGICAL LL1(ndomainsz,nlayer+1)
     109     
     110!--------------------------------------------------------------------------------
     111! 1. INITIALIZATION
     112!--------------------------------------------------------------------------------
     113! 100  CONTINUE  ! continue tag without source, maybe need delete in future
     114
     115      !* 1.1 COMPUTATIONAL CONSTANTS
    112116      kidia=1
    113           kfdia=klon
    114 
    115  110  CONTINUE
    116 C
    117       ILEVM1=KLEV-1
    118       ILEVM2=KLEV-2
    119       ILEVH =KLEV/3
    120 C
     117      kfdia=ngrid
     118! 110  CONTINUE  ! continue tag without source, maybe need delete in future
     119      ILEVM1=nlayer-1
     120      ILEVM2=nlayer-2
     121      ILEVH =nlayer/3   !!!! maybe not enough for Mars, need check later
    121122      ZCONS1=1./r
    122 cold  ZCONS2=G**2/CPD
    123123      ZCONS2=g**2/cpp
    124 cold  ZCONS3=1.5*API
    125124      ZCONS3=1.5*PI
    126 C
    127 C
    128 C     ------------------------------------------------------------------
    129 C
    130 C*         2.
    131 C                --------------
    132 C
    133  200  CONTINUE
    134 C
    135 C     ------------------------------------------------------------------
    136 C
    137 C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
    138 C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
    139 C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
    140 C
    141 C
    142 C
    143       DO 2001 JL=kidia,kfdia
    144       kknu(JL)    =klev
    145       kknu2(JL)   =klev
    146       kknub(JL)   =klev
    147       kknul(JL)   =klev
    148       pgamma(JL) =max(pgamma(jl),gtsec)
    149       ll1(jl,klev+1)=.false.
    150  2001 CONTINUE
    151 C
    152 C*      DEFINE TOP OF LOW LEVEL FLOW
    153 C       ----------------------------
    154       DO 2002 JK=KLEV,ilevh,-1
    155       DO 2003 JL=kidia,kfdia
    156       LO=(PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)).GE.GSIGCR
    157       IF(LO) THEN
    158         KKCRIT(JL)=JK
    159       ENDIF
    160       ZHCRIT(JL,JK)=4.*pvaror(JL)
    161       ZHGEO=PGEOM1(JL,JK)/g
    162       ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
    163       IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
    164         kknu(JL)=JK
    165       ENDIF
    166  2003 CONTINUE
    167  2002 CONTINUE
    168       DO 2004 JK=KLEV,ilevh,-1
    169       DO 2005 JL=kidia,kfdia
    170       ZHCRIT(JL,JK)=3.*pvaror(JL)
    171       ZHGEO=PGEOM1(JL,JK)/g
    172       ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
    173       IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
    174         kknu2(JL)=JK
    175       ENDIF
    176  2005 CONTINUE
    177  2004 CONTINUE
    178       DO 2006 JK=KLEV,ilevh,-1
    179       DO 2007 JL=kidia,kfdia
    180       ZHCRIT(JL,JK)=2.*pvaror(JL)
    181       ZHGEO=PGEOM1(JL,JK)/g
    182       ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
    183       IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
    184         kknub(JL)=JK
    185       ENDIF
    186  2007 CONTINUE
    187  2006 CONTINUE
    188       DO 2008 JK=KLEV,ilevh,-1
    189       DO 2009 JL=kidia,kfdia
    190       ZHCRIT(JL,JK)=pvaror(JL)
    191       ZHGEO=PGEOM1(JL,JK)/g
    192       ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
    193       IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
    194         kknul(JL)=JK
    195       ENDIF
    196  2009 CONTINUE
    197  2008 CONTINUE
    198 C
    199       do 2010 jl=kidia,kfdia 
    200       kknu(jl)=min(kknu(jl),nktopg)
    201       kknub(jl)=min(kknub(jl),nktopg)
    202       if(kknub(jl).eq.nktopg) kknul(jl)=klev
    203 C
    204 C     CHANGE IN HERE TO STOP KKNUL=KKNUB
    205 C
    206       if(kknul(jl).le.kknub(jl)) kknul(jl)=nktopg
    207  2010 continue     
    208 C
    209 
    210  210  CONTINUE
    211 C
    212 C
    213 CC*     INITIALIZE VARIOUS ARRAYS
    214 C
    215       DO 2107 JL=kidia,kfdia
    216       PRHO(JL,KLEV+1)  =0.0
    217       PSTAB(JL,KLEV+1) =0.0
    218       PSTAB(JL,1)      =0.0
    219       PRI(JL,KLEV+1)   =9999.0
    220       ppsi(JL,KLEV+1)  =0.0
    221       PRI(JL,1)        =0.0
    222       PVPH(JL,1)       =0.0
    223       PULOW(JL)        =0.0
    224       PVLOW(JL)        =0.0
    225       zulow(JL)        =0.0
    226       zvlow(JL)        =0.0
    227       KKCRITH(JL)      =KLEV
    228       KKenvH(JL)       =KLEV
    229       Kentp(JL)        =KLEV
    230       KCRIT(JL)        =1
    231       ncount(JL)       =0
    232       ll1(JL,klev+1)   =.false.
    233  2107 CONTINUE
    234 C
    235 C*     DEFINE LOW-LEVEL FLOW
    236 C      ---------------------
    237 C
    238       DO 223 JK=KLEV,2,-1
    239       DO 222 JL=kidia,kfdia
    240       IF(KTEST(JL).EQ.1) THEN
    241         ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
    242         PRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
    243         PSTAB(JL,JK)=2.*ZCONS2/(PTM1(JL,JK)+PTM1(JL,JK-1))*
    244      *  (1.-cpp*PRHO(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
    245         PSTAB(JL,JK)=MAX(PSTAB(JL,JK),GSSEC)
    246       ENDIF
    247   222 CONTINUE
    248   223 CONTINUE
    249 C
    250 C********************************************************************
    251 C
    252 C*     DEFINE blocked FLOW
    253 C      -------------------
    254       DO 2115 JK=klev,ilevh,-1
    255       DO 2116 JL=kidia,kfdia
    256       if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
    257         pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
    258         pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
    259       end if
    260  2116 CONTINUE
    261  2115 CONTINUE
    262       DO 2110 JL=kidia,kfdia
    263       pulow(JL)=pulow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
    264       pvlow(JL)=pvlow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
    265       ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
    266       PVPH(JL,KLEV+1)=ZNORM(JL)
    267  2110 CONTINUE
    268 C
    269 C*******  SETUP OROGRAPHY AXES AND DEFINE PLANE OF PROFILES  *******
    270 C
    271       DO 2112 JL=kidia,kfdia
    272       LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
    273       IF(LO) THEN
    274         ZU=PULOW(JL)+2.*GVSEC
    275       ELSE
    276         ZU=PULOW(JL)
    277       ENDIF
    278       Zphi=ATAN(PVLOW(JL)/ZU)
    279       ppsi(jl,klev+1)=ptheta(jl)*pi/180.-zphi
    280       zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
    281       zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
    282       pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
    283       pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
    284       pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
    285  2112 CONTINUE
    286 C
    287 C  ************ DEFINE FLOW IN PLANE OF LOWLEVEL STRESS *************
    288 C
    289       DO 213 JK=1,KLEV
    290       DO 212 JL=kidia,kfdia
    291       IF(KTEST(JL).EQ.1)  THEN
    292         ZVt1       =PULOW(JL)*PUM1(JL,JK)+PVLOW(JL)*PVM1(JL,JK)
    293         ZVt2       =-PvLOW(JL)*PUM1(JL,JK)+PuLOW(JL)*PVM1(JL,JK)
    294         ZVPF(JL,JK)=(zvt1*pd1(jl)+zvt2*pd2(JL))/(znorm(jl)*pdmod(jl))
    295       ENDIF
    296       PTAU(JL,JK)  =0.0
    297       Pzdep(JL,JK) =0.0
    298       Ppsi(JL,JK)  =0.0
    299       ll1(JL,JK)   =.FALSE.
    300   212 CONTINUE
    301   213 CONTINUE
    302       DO 215 JK=2,KLEV
    303       DO 214 JL=kidia,kfdia
    304       IF(KTEST(JL).EQ.1) THEN
    305         ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
    306         PVPH(JL,JK)=((PAPHM1(JL,JK)-PAPM1(JL,JK-1))*ZVPF(JL,JK)+
    307      *            (PAPM1(JL,JK)-PAPHM1(JL,JK))*ZVPF(JL,JK-1))
    308      *            /ZDP(JL,JK)
    309         IF(PVPH(JL,JK).LT.GVSEC) THEN
    310           PVPH(JL,JK)=GVSEC
    311           KCRIT(JL)=JK
     125
     126!------------------------------------------------------------------------------------------------------
     127! 2. Compute all the critical levels and the coeffecients of anisotropy
     128!-----------------------------------------------------------------------------------------------------
     129! 200  CONTINUE ! continue tag without source, maybe need delete in future
     130      ! 2.1 Define low level wind, project winds in plane of low level wind,
     131      ! determine sector in which to take the variance and set indicator for critical levels.
     132      DO JL=kidia,kfdia
     133           ! initialize all the height into surface (notice the layers have been inversed by preious rountines)
     134           IKNU(JL)    =nlayer
     135           IKNU2(JL)   =nlayer
     136           IKNUb(JL)   =nlayer
     137           IKNUl(JL)   =nlayer
     138           pgam(JL) =max(pgam(jl),gtsec)   ! gtsec is from yoegwd.h which is a common variable
     139           ll1(jl,nlayer+1)=.false.
     140      end DO
     141
     142      ! Define top of low level flow (since pressure, zonal and meridional wind have been inversed
     143      ! the process is to find the layer from surface (nlayer) to some levels ) by searching several
     144      ! altitude scope
     145     
     146      ! using 4 times sub-grid scale deviation as the threahold of the critical height
     147      DO JK=nlayer,ilevh,-1   ! ilevh=nlayer/3=16 
     148        DO JL=kidia,kfdia     ! jl=1:ngrid
     149           ! To found the layer of the "top of low level flow"
     150           LO=(pplev(JL,JK)/pplev(JL,nlayer+1)).GE.GSIGCR               
     151           IF(LO) THEN
     152             IKCRIT(JL)=JK
     153           ENDIF               
     154           ZHCRIT(JL,JK)=4.*pvar(JL) !
     155           ! use geopotetial denfination to get geoheight[in meters] of the layer
     156           ZHGEO=zgeom(JL,JK)/g     
     157           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))                 
     158           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
     159             IKNU(JL)=JK
     160           ENDIF
     161        ENDDO
     162      end DO
     163     
     164      ! using 3 times sub-grid scale deviation as the threahold of the critical height
     165      DO JK=nlayer,ilevh,-1
     166        DO JL=kidia,kfdia
     167           ZHCRIT(JL,JK)=3.*pvar(JL)
     168           ZHGEO=zgeom(JL,JK)/g
     169           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
     170           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
     171             IKNU2(JL)=JK
     172           ENDIF
     173        ENDDO
     174      end DO
     175     
     176      ! using 2 times sub-grid scale deviation as the threahold of the critical height
     177      DO JK=nlayer,ilevh,-1
     178        DO JL=kidia,kfdia
     179           ZHCRIT(JL,JK)=2.*pvar(JL)
     180           ZHGEO=zgeom(JL,JK)/g
     181           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
     182           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
     183             IKNUb(JL)=JK
     184           ENDIF
     185        ENDDO
     186      end DO
     187     
     188      ! using 1 times sub-grid scale deviation as the threahold of the critical height
     189      DO JK=nlayer,ilevh,-1
     190        DO JL=kidia,kfdia
     191           ZHCRIT(JL,JK)=pvar(JL)
     192           ZHGEO=zgeom(JL,JK)/g
     193           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
     194           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
     195             IKNUl(JL)=JK
     196           ENDIF
     197        ENDDO
     198      end DO     
     199      ! loop to relocate the critical height to make sure everything is okay if theses
     200      ! levels hit the model surface or top.
     201      do jl=kidia,kfdia 
     202         IKNU(jl)=min(IKNU(jl),nktopg)  ! nktopg is a common variable from yoegwd.h
     203         IKNUb(jl)=min(IKNUb(jl),nktopg)
     204         if(IKNUb(jl).eq.nktopg) IKNUl(jl)=nlayer
     205         ! Change in here to stop IKNUl=IKNUb
     206         if(IKNUl(jl).le.IKNUb(jl)) IKNUl(jl)=nktopg
     207      enddo     
     208
     209! 210  CONTINUE ! continue tag without source, maybe need delete in future
     210      ! Initialize various arrays for the following computes
     211      DO JL=kidia,kfdia
     212        ZRHO(JL,nlayer+1)  =0.0
     213        BV(JL,nlayer+1) =0.0
     214        BV(JL,1)      =0.0
     215        PRI(JL,nlayer+1)   =9999.0
     216        ZPSI(JL,nlayer+1)  =0.0
     217        PRI(JL,1)        =0.0
     218        ZVPH(JL,1)       =0.0
     219        PULOW(JL)        =0.0
     220        PVLOW(JL)        =0.0
     221        zulow(JL)        =0.0
     222        zvlow(JL)        =0.0
     223        IKCRITH(JL)      =nlayer     ! surface
     224        IKENVH(JL)       =nlayer     ! surface
     225        Kentp(JL)        =nlayer     ! surface
     226        ICRIT(JL)        =1          ! topmost layer
     227        ncount(JL)       =0
     228        ll1(JL,nlayer+1)   =.false.
     229      ENDDO
     230
     231      ! Define low-level flow Brunt–Väisälä frequency N^2, density ZRHO
     232      ! The incident flow passes over the mean orography is evaluated by averaging the wind,
     233      ! the Brunt–Väisälä frequency, and the fluid density between 1*pvar and 2*pvar over the
     234      ! model mean orography
     235      DO JK=nlayer,2,-1   ! from surface to topmost-1 layer
     236        DO JL=kidia,kfdia
     237           IF(ktest(JL).EQ.1) THEN ! if the map of the calling points is true
     238             ! calcalate density and BV
     239             ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1)  !dp>0
     240             ZRHO(JL,JK)=2.*pplev(JL,JK)*ZCONS1/(pt(JL,JK)+pt(JL,JK-1)) !rho=p/(r*T)
     241             !  Brunt–Väisälä frequency N^2. This equation for BV is illness since
     242             !  too many variables are used. Use N^2=g/T[1/(cpp*T)+dT/dz] to replace in the future
     243             BV(JL,JK)=2.*ZCONS2/(pt(JL,JK)+pt(JL,JK-1))*(1.-cpp*ZRHO(JL,JK)*(pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
     244             BV(JL,JK)=MAX(BV(JL,JK),GSSEC)
     245           ENDIF
     246        ENDDO
     247      end DO
     248     
     249      DO JK=nlayer,ilevh,-1
     250        DO JL=kidia,kfdia
     251           if(jk.ge.IKNUb(jl).and.jk.le.IKNUl(jl)) then ! IF the layer between 1*pvar and 2*pvar
     252           ! calculate the low level wind U_H
     253           ! pulow/pvlow at a speicfic location equals to sum of u*dp of all levels
     254           ! notice here dp is already a positive number
     255             pulow(JL)=pulow(JL)+pu(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
     256             pvlow(JL)=pvlow(JL)+pv(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
     257           end if
     258        ENDDO
     259      end DO
     260      ! averaging the wind
     261      DO JL=kidia,kfdia
     262         ! by divide dp [p differ between iknul and uknub level]
     263         pulow(JL)=pulow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
     264         pvlow(JL)=pvlow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
     265         ! average U to get background U?
     266         ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
     267         ZVPH(JL,nlayer+1)=ZNORM(JL)  ! The wind below the surface level (e.g., start of the 1/2 level)
     268      ENDDO     
     269     
     270      ! The gravity wave drag caused by the flow passes over an single elliptic mountain can be calculated
     271      ! by equation 17 and 18
     272      DO JL=kidia,kfdia
     273         LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
     274         IF(LO) THEN
     275           ZU=PULOW(JL)+2.*GVSEC
     276         ELSE
     277           ZU=PULOW(JL)
     278         ENDIF
     279         ! Here all physics for equation 17 and 18
     280         ! Direction of the incident flow
     281         Zphi=ATAN(PVLOW(JL)/ZU)
     282         ! The angle between the incident flow direction and the normal ridge direction pthe
     283         ZPSI(jl,nlayer+1)=pthe(jl)*pi/180.-zphi
     284         ! equation(17) parameter B and C
     285         zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2         
     286         zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2
     287         ! Bcos^2(psi)-Csin^2(psi)
     288         ZD1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ZPSI(jl,nlayer+1))**2)
     289         ! (B-C)sin(psi)cos(psi)
     290         ZD2(jl)=(zb(jl)-zc(jl))*sin(ZPSI(jl,nlayer+1))*cos(ZPSI(jl,nlayer+1))
     291         ! squre root of tao1 and tao2 without the constant see equation 17 or 18
     292         ZDMOD(jl)=sqrt(ZD1(jl)**2+ZD2(jl)**2)
     293      ENDDO
     294     
     295      ! Define blocked flow
     296      ! Setup orogrphy axes and define plane of profiles
     297      ! Define blocked flow in plane of the low level stress
     298      DO JK=1,nlayer
     299        DO JL=kidia,kfdia
     300           IF(ktest(JL).EQ.1)  THEN
     301             ZVt1       =PULOW(JL)*pu(JL,JK)+PVLOW(JL)*pv(JL,JK)
     302             ZVt2       =-PvLOW(JL)*pu(JL,JK)+PuLOW(JL)*pv(JL,JK)
     303             ! zvpf is a normalized variable
     304             ZVPF(JL,JK)=(zvt1*ZD1(jl)+zvt2*ZD2(JL))/(znorm(jl)*ZDMOD(jl))
     305           ENDIF
     306           ZTAU(JL,JK)  =0.0
     307           ZZDEP(JL,JK) =0.0
     308           ZPSI(JL,JK)  =0.0
     309           ll1(JL,JK)   =.FALSE.
     310        ENDDO
     311      end DO
     312 
     313      DO JK=2,nlayer
     314        DO JL=kidia,kfdia
     315           IF(ktest(JL).EQ.1) THEN
     316             ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1)  ! dp
     317             ! zvph is the U_H in equation 17 e.g. low level wind speed
     318             ZVPH(JL,JK)=((pplev(JL,JK)-pplay(JL,JK-1))*ZVPF(JL,JK)+ &
     319             (pplay(JL,JK)-pplev(JL,JK))*ZVPF(JL,JK-1))/ZDP(JL,JK)
     320             IF(ZVPH(JL,JK).LT.GVSEC) THEN
     321               ZVPH(JL,JK)=GVSEC
     322               ICRIT(JL)=JK
     323             ENDIF
     324           endIF
     325        ENDDO
     326      end DO
     327
     328      ! 2.2 Brunt-vaisala frequency and density at half levels
     329 220  CONTINUE ! continue tag without source, maybe need delete in future
     330 
     331      DO JK=ilevh,nlayer
     332        DO JL=kidia,kfdia
     333           IF(ktest(JL).EQ.1) THEN
     334              IF(jk.ge.(IKNUb(jl)+1).and.jk.le.IKNUl(jl)) THEN
     335                 ZST=ZCONS2/pt(JL,JK)*(1.-cpp*ZRHO(JL,JK)*     &
     336                 (pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
     337                 BV(JL,nlayer+1)=BV(JL,nlayer+1)+ZST*ZDP(JL,JK)
     338                 BV(JL,nlayer+1)=MAX(BV(JL,nlayer+1),GSSEC)
     339                 ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)+pplev(JL,JK)*2.*ZDP(JL,JK) &
     340                 *ZCONS1/(pt(JL,JK)+pt(JL,JK-1))
     341              ENDIF
     342           endIF
     343        ENDDO
     344      end DO
     345
     346      DO JL=kidia,kfdia
     347!*****************************************************************************
     348!     Okay. There is a possible problem here. If IKNUl=IKNUb then division by zero
     349!     occurs. I have put a fix in here but will ask Francois lott about it in Paris.     
     350!     Also if this is the case BV and ZRHO are not defined at nlayer+1 so I have
     351!     added the else.
     352!     by: MAT COLLINS 30.1.96
     353!*****************************************************************************
     354        IF (IKNUL(JL).NE.IKNUB(JL)) THEN
     355           BV(JL,nlayer+1)=BV(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
     356           ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
     357        ELSE
     358           WRITE(*,*) 'OROSETUP: IKNUB=IKNUL= ',IKNUB(JL),' AT JL= ',JL
     359           BV(JL,nlayer+1)=BV(JL,nlayer)
     360           ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer)
    312361        ENDIF
    313       ENDIF
    314   214 CONTINUE
    315   215 CONTINUE
    316 C
    317 C
    318 C*         2.2     BRUNT-VAISALA FREQUENCY AND DENSITY AT HALF LEVELS.
    319 C
    320   220 CONTINUE
    321 C
    322       DO 2211 JK=ilevh,KLEV
    323       DO 221 JL=kidia,kfdia
    324       IF(KTEST(JL).EQ.1) THEN
    325       IF(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) THEN
    326            ZST=ZCONS2/PTM1(JL,JK)*(1.-cpp*PRHO(JL,JK)*
    327      *                   (PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
    328            PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)+ZST*ZDP(JL,JK)
    329            PSTAB(JL,KLEV+1)=MAX(PSTAB(JL,KLEV+1),GSSEC)
    330            PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)+PAPHM1(JL,JK)*2.*ZDP(JL,JK)
    331      *                   *ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
    332       ENDIF
    333       ENDIF
    334   221 CONTINUE
    335  2211 CONTINUE
    336 C
    337       DO 2212 JL=kidia,kfdia
    338 C*****************************************************************************
    339 C
    340 C     O.K. THERE IS A POSSIBLE PROBLEM HERE. IF KKNUL=KKNUB THEN
    341 C     DIVISION BY ZERO OCCURS. I HAVE PUT A FIX IN HERE BUT WILL ASK FRANCOIS
    342 C     LOTT ABOUT IT IN PARIS.
    343 C
    344 C     MAT COLLINS 30.1.96
    345 C
    346 C     ALSO IF THIS IS THE CASE PSTAB AND PRHO ARE NOT DEFINED AT KLEV+1
    347 C     SO I HAVE ADDED THE ELSE
    348 C
    349 C*****************************************************************************
    350       IF (KKNUL(JL).NE.KKNUB(JL)) THEN
    351         PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
    352      *                                          -PAPM1(JL,kknub(jl)))
    353         PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
    354      *                                          -PAPM1(JL,kknub(jl)))
    355       ELSE
    356       WRITE(*,*) 'OROSETUP: KKNUB=KKNUL= ',KKNUB(JL),' AT JL= ',JL
    357         PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV)
    358         PRHO(JL,KLEV+1)=PRHO(JL,KLEV)
    359       ENDIF
    360         ZVAR=PVARor(JL)
    361  2212 CONTINUE
    362 C
    363 C*         2.3     MEAN FLOW RICHARDSON NUMBER.
    364 C*                 AND CRITICAL HEIGHT FOR FROUDE LAYER
    365 C
    366   230 CONTINUE
    367 C
    368       DO 232 JK=2,KLEV
    369       DO 231 JL=kidia,kfdia
    370       IF(KTEST(JL).EQ.1) THEN
    371         ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
    372         PRI(JL,JK)=PSTAB(JL,JK)*(ZDP(JL,JK)
    373      *          /(g*PRHO(JL,JK)*ZDWIND))**2
    374         PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
    375       ENDIF
    376   231 CONTINUE
    377   232 CONTINUE
    378 C
    379 C
    380 C*      DEFINE TOP OF 'envelope' layer
    381 C       ----------------------------
    382 
    383       DO 233 JL=kidia,kfdia
    384       pnu (jl)=0.0
    385       znum(jl)=0.0
    386  233  CONTINUE
    387      
    388       DO 234 JK=2,KLEV-1
    389       DO 234 JL=kidia,kfdia
    390      
    391       IF(KTEST(JL).EQ.1) THEN
    392        
    393       IF (JK.GE.KKNU2(JL)) THEN
    394          
    395             ZNUM(JL)=PNU(JL)
    396             ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
    397      *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
    398             ZWIND=max(sqrt(zwind**2),gvsec)
    399             ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
    400             ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
    401             ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
    402             ZRHOM=PRHO(JL,JK  )
    403             ZRHOP=PRHO(JL,JK+1)
    404             PNU(JL) = PNU(JL) + (ZDELP/g)*
    405      *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
    406             IF((ZNUM(JL).LE.GFRCRIT).AND.(PNU(JL).GT.GFRCRIT)
    407      *                          .AND.(KKENVH(JL).EQ.KLEV))
    408      *      KKENVH(JL)=JK
    409      
    410       ENDIF   
    411 
    412       ENDIF
    413      
    414  234  CONTINUE
    415 
    416 C  CALCULATION OF A DYNAMICAL MIXING HEIGHT FOR THE BREAKING
    417 C  OF GRAVITY WAVES:
    418 
    419              
    420       DO 235 JL=kidia,kfdia
    421       znup(jl)=0.0
    422       znum(jl)=0.0
    423  235  CONTINUE
    424 
    425       DO 236 JK=KLEV-1,2,-1
    426       DO 236 JL=kidia,kfdia
    427      
    428       IF(KTEST(JL).EQ.1) THEN
    429        
    430       IF (JK.LT.KKENVH(JL)) THEN
    431 
    432             ZNUM(JL)=ZNUP(JL)
    433             ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
    434      *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
    435             ZWIND=max(sqrt(zwind**2),gvsec)
    436             ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
    437             ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
    438             ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
    439             ZRHOM=PRHO(JL,JK  )
    440             ZRHOP=PRHO(JL,JK+1)
    441             ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*
    442      *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
    443             IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5)
    444      *                          .AND.(KKCRITH(JL).EQ.KLEV))
    445      *      KKCRITH(JL)=JK
    446      
    447       ENDIF
    448      
    449       ENDIF
    450      
    451  236  CONTINUE
    452  
    453       DO 237 JL=KIDIA,KFDIA
    454       KKCRITH(JL)=MIN0(KKCRITH(JL),KKNU(JL))
    455  237  CONTINUE
    456 c
    457 c     directional info for flow blocking *************************
    458 c
    459       do 251 jk=ilevh,klev   
    460       DO 252 JL=kidia,kfdia
    461       IF(jk.ge.kkenvh(jl)) THEN
    462       LO=(PUm1(JL,jk).LT.GVSEC).AND.(PUm1(JL,jk).GE.-GVSEC)
    463       IF(LO) THEN
    464         ZU=PUm1(JL,jk)+2.*GVSEC
    465       ELSE
    466         ZU=PUm1(JL,jk)
    467       ENDIF
    468        Zphi=ATAN(PVm1(JL,jk)/ZU)
    469        ppsi(jl,jk)=ptheta(jl)*pi/180.-zphi
    470       end if
    471  252  continue
    472  251  CONTINUE
    473 c      forms the vertical 'leakiness' **************************
    474 
    475       alpha=3.
    476      
    477       DO 254  JK=ilevh,klev
    478       DO 253  JL=kidia,kfdia
    479       IF(jk.ge.kkenvh(jl)) THEN
    480         zggeenv=AMAX1(1.,
    481      *          (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)     
    482         zggeom1=AMAX1(pgeom1(jl,jk),1.)
    483         zgvar=amax1(pvaror(jl)*g,1.)     
    484         pzdep(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))     
    485       end if
    486  253  CONTINUE
    487  254  CONTINUE
    488 
    489  260  CONTINUE
    490 
    491       RETURN
    492       END
     362        ZVAR=pvar(JL)
     363      ENDDO !JL=kidia,kfdia
     364
     365      ! 2.3 Mean flow richardson number and critical height for proude layer   
     366! 230  CONTINUE ! continue tag without source, maybe need delete in future
     367
     368      DO JK=2,nlayer
     369        DO JL=kidia,kfdia
     370           IF(ktest(JL).EQ.1) THEN
     371             ! du
     372             ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
     373             ! Mean flow Richardson number Ri=g/rho[drho/dz / (du/dz)^2]
     374             ! Here dp maybe dp^2 ? Need ask Francios lott later
     375             PRI(JL,JK)=BV(JL,JK)*(ZDP(JL,JK)/(g*ZRHO(JL,JK)*ZDWIND))**2
     376             PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
     377           ENDIF
     378        ENDDO
     379      end DO
     380
     381      !*    Define top of 'envelope' layer 
     382      DO JL=kidia,kfdia
     383         ZNU (jl)=0.0
     384         znum(jl)=0.0
     385      ENDDO
     386     
     387      DO JK=2,nlayer-1
     388        DO JL=kidia,kfdia     
     389           IF(ktest(JL).EQ.1) THEN       
     390             IF (JK.GE.IKNU2(JL)) THEN  ! level lower than 3*par
     391             ! all codes here is to calculate equation 9       
     392                ZNUM(JL)=ZNU(JL)
     393                ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
     394                ZWIND=max(sqrt(zwind**2),gvsec)
     395                ZDELP=pplev(JL,JK+1)-pplev(JL,JK)   ! dp
     396                ZSTABM=SQRT(MAX(BV(JL,JK  ),GSSEC))
     397                ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
     398                ZRHOM=ZRHO(JL,JK  )
     399                ZRHOP=ZRHO(JL,JK+1)
     400                ! Equation 9. znu is a critical value to find the blocking layer
     401                ZNU(JL) = ZNU(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND 
     402                ! Found the moutain top   
     403                IF((ZNUM(JL).LE.GFRCRIT).AND.(ZNU(JL).GT.GFRCRIT).AND.(IKENVH(JL).EQ.nlayer)) THEN
     404                  IKENVH(JL)=JK
     405                ENDIF     
     406             ENDIF ! (JK.GE.IKNU2(JL)) 
     407           ENDIF !(ktest(JL).EQ.1)   
     408        ENDDO
     409       endDO
     410
     411      !  Calculation of a dynamical mixing height for the breaking of gravity waves           
     412      DO JL=kidia,kfdia
     413         znup(jl)=0.0
     414         znum(jl)=0.0
     415      ENDDO
     416
     417      DO JK=nlayer-1,2,-1
     418        DO JL=kidia,kfdia     
     419          IF(ktest(JL).EQ.1) THEN       
     420            IF (JK.LT.IKENVH(JL)) THEN
     421                ZNUM(JL)=ZNUP(JL)
     422                ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
     423                ZWIND=max(sqrt(zwind**2),gvsec)
     424                ZDELP=pplev(JL,JK+1)-pplev(JL,JK)
     425                ZSTABM=SQRT(MAX(BV(JL,JK  ),GSSEC))
     426                ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
     427                ZRHOM=ZRHO(JL,JK  )
     428                ZRHOP=ZRHO(JL,JK+1)
     429                ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND 
     430                ! dynamical mixing height for the breaking of gravity waves   
     431                IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5).AND.(IKCRITH(JL).EQ.nlayer)) THEN
     432                   IKCRITH(JL)=JK
     433                ENDIF     
     434            ENDIF  ! (JK.LT.IKENVH(JL))   
     435          ENDIF    ! (ktest(JL).EQ.1)   
     436        ENDDO
     437      end DO
     438     
     439      DO JL=KIDIA,KFDIA
     440         IKCRITH(JL)=MIN0(IKCRITH(JL),IKNU(JL))
     441      ENDDO
     442
     443      ! directional info for flow blocking *************************
     444      DO jk=ilevh,nlayer   
     445        DO JL=kidia,kfdia
     446           IF(jk.ge.IKENVH(jl)) THEN
     447             LO=(pu(JL,jk).LT.GVSEC).AND.(pu(JL,jk).GE.-GVSEC)
     448             IF(LO) THEN
     449               ZU=pu(JL,jk)+2.*GVSEC
     450             ELSE
     451               ZU=pu(JL,jk)
     452             ENDIF
     453             Zphi=ATAN(pv(JL,jk)/ZU)
     454             ZPSI(jl,jk)=pthe(jl)*pi/180.-zphi
     455           end IF
     456        ENDDO
     457      end DO
     458     
     459      ! forms the vertical 'leakiness' ************************** 
     460      DO JK=ilevh,nlayer
     461        DO JL=kidia,kfdia
     462           IF(jk.ge.IKENVH(jl)) THEN
     463             zggeenv=AMAX1(1.,(zgeom(jl,IKENVH(jl))+zgeom(jl,IKENVH(jl)-1))/2.)     
     464             zggeom1=AMAX1(zgeom(jl,jk),1.)
     465             zgvar=amax1(pvar(jl)*g,1.)     
     466             ZZDEP(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))     
     467          endIF
     468        ENDDO
     469      end DO
     470
     471! 260  CONTINUE  ! continue tag without source, maybe need delete in future
     472
     473  RETURN
     474END
Note: See TracChangeset for help on using the changeset viewer.