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

File:
1 moved

Legend:

Unmodified
Added
Removed

  • 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 
Note: See TracChangeset for help on using the changeset viewer.