Ignore:
Timestamp:
Jun 26, 2014, 6:07:05 PM (11 years ago)
Author:
emillour
Message:

Common dynamics:
Some updates to keep up with LMDZ5 Earth model evolution
(up to LMDZ5 rev 2070). See file "DOC/chantiers/commit_importants.log"
for detailed list of changes.
Note that the updates of exner* routines change (as expected) results
at numerical roundoff level.
EM

Location:
trunk/LMDZ.COMMON/libf/dyn3d_common
Files:
3 deleted
4 edited
2 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90

    r1300 r1302  
    11! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
    22
    3 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
     3SUBROUTINE disvert()
    44
    55  ! Auteur : P. Le Van
    66
     7#ifdef CPP_IOIPSL
     8  use ioipsl, only: getin
     9#else
     10  USE ioipsl_getincom, only: getin
     11#endif
    712  use new_unit_m, only: new_unit
    8   use ioipsl, only: getin
    913  use assert_m, only: assert
    1014
    … …  
    1317  include "dimensions.h"
    1418  include "paramet.h"
     19  include "comvert.h"
     20  include "comconst.h"
    1521  include "iniprint.h"
    1622  include "logic.h"
    1723
    18   ! s = sigma ** kappa : coordonnee verticale
    19   ! dsig(l) : epaisseur de la couche l ds la coord. s
    20   ! sig(l) : sigma a l'interface des couches l et l-1
    21   ! ds(l) : distance entre les couches l et l-1 en coord.s
    22 
    23   real,intent(in) :: pa, preff
    24   real,intent(out) :: ap(llmp1) ! in Pa
    25   real, intent(out):: bp(llmp1)
    26   real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
    27   real,intent(out) :: presnivs(llm)
    28   real,intent(out) :: scaleheight
    29 
     24!-------------------------------------------------------------------------------
     25! Purpose: Vertical distribution functions for LMDZ.
     26!          Triggered by the levels number llm.
     27!-------------------------------------------------------------------------------
     28! Read    in "comvert.h":
     29! pa                         !--- PURE PRESSURE COORDINATE FOR P<pa (in Pascals)
     30! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
     31! Written in "comvert.h":
     32! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
     33! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
     34! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
     35! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
     36! scaleheight                !--- VERTICAL SCALE HEIGHT            (Earth: 8kms)
     37! nivsig(llm+1)              !--- SIGMA INDEX OF EACH LAYER INTERFACE
     38! nivsigs(llm)               !--- SIGMA INDEX OF EACH MID-LAYER
     39!-------------------------------------------------------------------------------
     40! Local variables:
    3041  REAL sig(llm+1), dsig(llm)
    31   real zk, zkm1, dzk1, dzk2, k0, k1
     42  REAL sig0(llm+1), zz(llm+1)
     43  REAL zk, zkm1, dzk1, dzk2, z, k0, k1
    3244
    3345  INTEGER l, unit
    3446  REAL dsigmin
     47  REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
     48
    3549  REAL alpha, beta, deltaz
    3650  REAL x
    3751  character(len=*),parameter :: modname="disvert"
    3852
    39   character(len=6):: vert_sampling
     53  character(len=24):: vert_sampling
    4054  ! (allowed values are "param", "tropo", "strato" and "read")
    4155
    4256  !-----------------------------------------------------------------------
    4357
    44   print *, "Call sequence information: disvert"
     58  WRITE(lunout,*) TRIM(modname)//" starts"
    4559
    4660  ! default scaleheight is 8km for earth
    … …  
    4963  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
    5064  call getin('vert_sampling', vert_sampling)
    51   print *, 'vert_sampling = ' // vert_sampling
     65  WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
     66  if (llm==39 .and. vert_sampling=="strato") then
     67     dsigmin=0.3 ! Vieille option par défaut pour CMIP5
     68  else
     69     dsigmin=1.
     70  endif
     71  call getin('dsigmin', dsigmin)
     72  WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
     73
    5274
    5375  select case (vert_sampling)
    … …  
    86108
    87109     ap = pa * (sig - bp)
    88   case("tropo")
     110  case("sigma")
    89111     DO l = 1, llm
    90112        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    91         dsig(l) = 1.0 + 7.0 * SIN(x)**2
     113        dsig(l) = dsigmin + 7.0 * SIN(x)**2
    92114     ENDDO
    93115     dsig = dsig / sum(dsig)
    … …  
    98120
    99121     bp(1)=1.
     122     bp(2: llm) = sig(2:llm)
     123     bp(llmp1) = 0.
     124     ap(:)=0.
     125  case("tropo")
     126     DO l = 1, llm
     127        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     128        dsig(l) = dsigmin + 7.0 * SIN(x)**2
     129     ENDDO
     130     dsig = dsig / sum(dsig)
     131     sig(llm+1) = 0.
     132     DO l = llm, 1, -1
     133        sig(l) = sig(l+1) + dsig(l)
     134     ENDDO
     135
     136     bp(1)=1.
    100137     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
    101138     bp(llmp1) = 0.
    … …  
    104141     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
    105142  case("strato")
    106      if (llm==39) then
    107         dsigmin=0.3
    108      else if (llm==50) then
    109         dsigmin=1.
    110      else
    111         write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
    112         dsigmin=1.
    113      endif
    114      WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
    115 
    116143     DO l = 1, llm
    117144        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    … …  
    131158     ap(1)=0.
    132159     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
     160  case("strato_correct")
     161!==================================================================
     162! Fredho 2014/05/18, Saint-Louis du Senegal
     163! Cette version de la discretisation strato est corrige au niveau
     164! du passage des sig aux ap, bp
     165! la version precedente donne un coude dans l'epaisseur des couches
     166! vers la tropopause
     167!==================================================================
     168
     169
     170     DO l = 1, llm
     171        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     172        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     173             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     174     ENDDO
     175     dsig = dsig / sum(dsig)
     176     sig0(llm+1) = 0.
     177     DO l = llm, 1, -1
     178        sig0(l) = sig0(l+1) + dsig(l)
     179     ENDDO
     180     sig=racinesig(sig0)
     181
     182     bp(1)=1.
     183     bp(2:llm)=EXP(1.-1./sig(2: llm)**2)
     184     bp(llmp1)=0.
     185
     186     ap(1)=0.
     187     ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
     188     ap(llm+1)=0.
     189
     190  CASE("strato_custom0")
     191!=======================================================
     192! Version Transitoire
     193    ! custumize strato distribution with specific alpha & beta values and function
     194    ! depending on llm (experimental and temporary)!
     195    SELECT CASE (llm)
     196      CASE(55)
     197        alpha=0.45
     198        beta=4.0
     199      CASE(63)
     200        alpha=0.45
     201        beta=5.0
     202      CASE(71)
     203        alpha=3.05
     204        beta=65.
     205      CASE(79)
     206        alpha=3.20
     207        ! alpha=2.05 ! FLOTT 79 (PLANTE)
     208        beta=70.
     209    END SELECT
     210    ! Or used values provided by user in def file:
     211    CALL getin("strato_alpha",alpha)
     212    CALL getin("strato_beta",beta)
     213   
     214    ! Build geometrical distribution
     215    scaleheight=7.
     216    zz(1)=0.
     217    IF (llm==55.OR.llm==63) THEN
     218      DO l=1,llm
     219        z=zz(l)/scaleheight
     220        zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5))     &
     221                            +5.0*EXP((l-llm)/beta)
     222      ENDDO
     223    ELSEIF (llm==71) THEN !.OR.llm==79) THEN      ! FLOTT 79 (PLANTE)
     224      DO l=1,llm
     225        z=zz(l)
     226        zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
     227      ENDDO
     228    ELSEIF (llm==79) THEN
     229      DO l=1,llm
     230        z=zz(l)
     231        zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.))     &
     232                            +0.03*TANH(z/.25)
     233      ENDDO
     234    ENDIF ! of IF (llm==55.OR.llm==63) ...
     235   
     236
     237    ! Build sigma distribution
     238    sig0=EXP(-zz(:)/scaleheight)
     239    sig0(llm+1)=0.
     240!    sig=ridders(sig0)
     241    sig=racinesig(sig0)
     242   
     243    ! Compute ap() and bp()
     244    bp(1)=1.
     245    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
     246    bp(llm+1)=0.
     247    ap=pa*(sig-bp)
     248
     249  CASE("strato_custom")
     250!===================================================================
     251! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
     252! 2014/05
     253! custumize strato distribution
     254! Al the parameter are given in km assuming a given scalehigh
     255    vert_scale=7.     ! scale hight
     256    vert_dzmin=0.02   ! width of first layer
     257    vert_dzlow=1.     ! dz in the low atmosphere
     258    vert_z0low=8.     ! height at which resolution recches dzlow
     259    vert_dzmid=3.     ! dz in the mid atmsophere
     260    vert_z0mid=70.    ! height at which resolution recches dzmid
     261    vert_h_mid=20.    ! width of the transition
     262    vert_dzhig=11.    ! dz in the high atmsophere
     263    vert_z0hig=80.    ! height at which resolution recches dz
     264    vert_h_hig=20.    ! width of the transition
     265!===================================================================
     266
     267    call getin('vert_scale',vert_scale)
     268    call getin('vert_dzmin',vert_dzmin)
     269    call getin('vert_dzlow',vert_dzlow)
     270    call getin('vert_z0low',vert_z0low)
     271    CALL getin('vert_dzmid',vert_dzmid)
     272    CALL getin('vert_z0mid',vert_z0mid)
     273    call getin('vert_h_mid',vert_h_mid)
     274    call getin('vert_dzhig',vert_dzhig)
     275    call getin('vert_z0hig',vert_z0hig)
     276    call getin('vert_h_hig',vert_h_hig)
     277
     278    scaleheight=vert_scale ! for consistency with further computations
     279    ! Build geometrical distribution
     280    zz(1)=0.
     281    DO l=1,llm
     282       z=zz(l)
     283       zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+                &
     284&      (vert_dzmid-vert_dzlow)* &
     285&           (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
     286&     +(vert_dzhig-vert_dzmid-vert_dzlow)*                                  &
     287&           (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
     288    ENDDO
     289
     290
     291!===================================================================
     292! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
     293! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
     294! sig0 is p/p0
     295! sig is an intermediate distribution introduce to estimate bp
     296! 1.   sig0=exp(-z/H)
     297! 2.   inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
     298! 3.   bp=exp(1-1/sig**2)
     299! 4.   ap deduced from  the combination of 2 and 3 so that sig0=ap/p0+bp
     300!===================================================================
     301
     302    sig0=EXP(-zz(:)/vert_scale)
     303    sig0(llm+1)=0.
     304    sig=racinesig(sig0)
     305    bp(1)=1.
     306    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
     307    bp(llm+1)=0.
     308    ap=pa*(sig-bp)
     309
    133310  case("read")
    134311     ! Read "ap" and "bp". First line is skipped (title line). "ap"
    … …  
    178355  write(lunout, *) presnivs
    179356
     357CONTAINS
     358
     359!-------------------------------------------------------------------------------
     360!
     361FUNCTION ridders(sig) RESULT(sg)
     362!
     363!-------------------------------------------------------------------------------
     364  IMPLICIT NONE
     365!-------------------------------------------------------------------------------
     366! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
     367! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
     368! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
     369!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
     370!-------------------------------------------------------------------------------
     371! Arguments:
     372  REAL, INTENT(IN)  :: sig(:)
     373  REAL              :: sg(SIZE(sig))
     374!-------------------------------------------------------------------------------
     375! Local variables:
     376  INTEGER :: it, ns, maxit
     377  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
     378!-------------------------------------------------------------------------------
     379  ns=SIZE(sig); maxit=9999
     380  c1=Pa/Preff; c2=1.-c1
     381  DO l=1,ns
     382    xx=HUGE(1.)
     383    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
     384    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
     385    DO it=1,maxit
     386      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
     387      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
     388      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
     389      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
     390      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
     391      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
     392      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
     393      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
     394      END IF
     395      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m           
     396    END DO
     397    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
     398     &e for level ',l
     399    sg(l)=xx
     400  END DO
     401  sg(1)=1.; sg(ns)=0.
     402
     403END FUNCTION ridders
     404
     405FUNCTION racinesig(sig) RESULT(sg)
     406!
     407!-------------------------------------------------------------------------------
     408  IMPLICIT NONE
     409!-------------------------------------------------------------------------------
     410! Fredho 2014/05/18
     411! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
     412! Notes:   Uses Newton Raphson search
     413!-------------------------------------------------------------------------------
     414! Arguments:
     415  REAL, INTENT(IN)  :: sig(:)
     416  REAL              :: sg(SIZE(sig))
     417!-------------------------------------------------------------------------------
     418! Local variables:
     419  INTEGER :: it, ns, maxit
     420  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
     421!-------------------------------------------------------------------------------
     422  ns=SIZE(sig); maxit=100
     423  c1=Pa/Preff; c2=1.-c1
     424  DO l=2,ns-1
     425    sg(l)=sig(l)
     426    DO it=1,maxit
     427       f1=exp(1-1./sg(l)**2)*(1.-c1)
     428       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
     429    ENDDO
     430!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
     431  ENDDO
     432  sg(1)=1.; sg(ns)=0.
     433
     434END FUNCTION racinesig
     435
     436
     437
     438
    180439END SUBROUTINE disvert
     440
     441!-------------------------------------------------------------------------------
     442
     443FUNCTION distrib(x,c1,c2,x0) RESULT(res)
     444!
     445!-------------------------------------------------------------------------------
     446! Arguments:
     447  REAL, INTENT(IN) :: x, c1, c2, x0
     448  REAL             :: res
     449!-------------------------------------------------------------------------------
     450  res=c1*x+c2*EXP(1-1/(x**2))-x0
     451
     452END FUNCTION distrib
     453
     454
     455

  • trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90

    r1300 r1302  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  P.Le Van  , Fr. Hourdin  .
    7 c    ..........
    8 c
    9 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    10 c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    11 c
    12 c   ************************************************************************
    13 c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
    14 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    15 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    16 c   ************************************************************************
    17 c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    18 c    la pression et la fonction d'Exner  au  sol  .
    19 c
    20 c                                 -------- z                                   
    21 c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
    22 c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
    23 c    ( voir note de Fr.Hourdin )  ,
    24 c
    25 c    on determine successivement , du haut vers le bas des couches, les
    26 c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
    27 c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
    28 c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
    29 c
    30 c
    31       IMPLICIT NONE
    32 c
    33 #include "dimensions.h"
    34 #include "paramet.h"
    35 #include "comconst.h"
    36 #include "comgeom.h"
    37 #include "comvert.h"
    38 #include "serre.h"
     1module exner_hyb_m
    392
    40       INTEGER  ngrid
    41       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    42       REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
     3  IMPLICIT NONE
    434
    44 c    .... variables locales   ...
     5contains
    456
    46       INTEGER l, ij
    47       REAL unpl2k,dellta
     7  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
    488
    49       REAL ppn(iim),pps(iim)
    50       REAL xpn, xps
    51       REAL SSUM
    52 c
    53       logical,save :: firstcall=.true.
    54       character(len=*),parameter :: modname="exner_hyb"
    55      
    56       ! Sanity check
    57       if (firstcall) then
    58         ! sanity checks for Shallow Water case (1 vertical layer)
    59         if (llm.eq.1) then
     9    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
     10    !    ..........
     11    !
     12    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     13    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     14    !
     15    !   ************************************************************************
     16    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     17    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     18    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     19    !   ************************************************************************
     20    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     21    !    la pression et la fonction d'Exner  au  sol  .
     22    !
     23    !                                 -------- z
     24    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
     25    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
     26    !    ( voir note de Fr.Hourdin )  ,
     27    !
     28    !    on determine successivement , du haut vers le bas des couches, les
     29    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
     30    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
     31    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
     32    !
     33    !
     34    !
     35    include "dimensions.h"
     36    include "paramet.h"
     37    include "comconst.h"
     38    include "comgeom.h"
     39    include "comvert.h"
     40    include "serre.h"
     41
     42    INTEGER  ngrid
     43    REAL p(ngrid,llmp1),pk(ngrid,llm)
     44    real, optional:: pkf(ngrid,llm)
     45    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
     46
     47    !    .... variables locales   ...
     48
     49    INTEGER l, ij
     50    REAL unpl2k,dellta
     51
     52    logical,save :: firstcall=.true.
     53    character(len=*),parameter :: modname="exner_hyb"
     54
     55    ! Sanity check
     56    if (firstcall) then
     57       ! sanity checks for Shallow Water case (1 vertical layer)
     58       if (llm.eq.1) then
    6059          if (kappa.ne.1) then
    61             call abort_gcm(modname,
    62      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     60             call abort_gcm(modname, &
     61                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6362          endif
    6463          if (cpp.ne.r) then
    65             call abort_gcm(modname,
    66      &      "cpp!=r , but running in Shallow Water mode!!",42)
     64             call abort_gcm(modname, &
     65                  "cpp!=r , but running in Shallow Water mode!!",42)
    6766          endif
    68         endif ! of if (llm.eq.1)
     67       endif ! of if (llm.eq.1)
    6968
    70         firstcall=.false.
    71       endif ! of if (firstcall)
     69       firstcall=.false.
     70    endif ! of if (firstcall)
    7271
    73       if (llm.eq.1) then
    74        
    75         ! Compute pks(:),pk(:),pkf(:)
    76        
    77         DO   ij  = 1, ngrid
    78           pks(ij) = (cpp/preff) * ps(ij)
     72    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     73    if (llm.eq.1) then
     74
     75       ! Compute pks(:),pk(:),pkf(:)
     76
     77       DO   ij  = 1, ngrid
     78          pks(ij) = (cpp/preff) * ps(ij)
    7979          pk(ij,1) = .5*pks(ij)
    80         ENDDO
    81        
    82         CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    83         CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    84        
    85         ! our work is done, exit routine
    86         return
     80       ENDDO
    8781
    88       endif ! of if (llm.eq.1)
     82       if (present(pkf)) then
     83          pkf = pk
     84          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     85       end if
    8986
    90 !!!! General case:
    91      
    92       unpl2k    = 1.+ 2.* kappa
    93 c
    94       DO   ij  = 1, ngrid
    95         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    96       ENDDO
     87       ! our work is done, exit routine
     88       return
     89    endif ! of if (llm.eq.1)
    9790
    98       DO  ij   = 1, iim
    99         ppn(ij) = aire(   ij   ) * pks(  ij     )
    100         pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    101       ENDDO
    102       xpn      = SSUM(iim,ppn,1) /apoln
    103       xps      = SSUM(iim,pps,1) /apols
     91    ! General case:
    10492
    105       DO ij   = 1, iip1
    106         pks(   ij     )  =  xpn
    107         pks( ij+ip1jm )  =  xps
    108       ENDDO
    109 c
    110 c
    111 c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
    112 c
    113       DO     ij      = 1, ngrid
     93    unpl2k    = 1.+ 2.* kappa
     94
     95    !     -------------
     96    !     Calcul de pks
     97    !     -------------
     98
     99    DO   ij  = 1, ngrid
     100       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     101    ENDDO
     102
     103    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
     104    !
     105    DO     ij      = 1, ngrid
    114106       alpha(ij,llm) = 0.
    115107       beta (ij,llm) = 1./ unpl2k
    116       ENDDO
    117 c
    118 c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
    119 c
    120       DO l = llm -1 , 2 , -1
    121 c
    122         DO ij = 1, ngrid
    123         dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
    124         alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
    125         beta (ij,l)  =   p(ij,l  ) / dellta   
    126         ENDDO
    127 c
    128       ENDDO
    129 c
    130 c  ***********************************************************************
    131 c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
    132 c
    133       DO   ij   = 1, ngrid
    134        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
    135      *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
    136       ENDDO
    137 c
    138 c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
    139 c
    140       DO l = 2, llm
    141         DO   ij   = 1, ngrid
    142          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
    143         ENDDO
    144       ENDDO
    145 c
    146 c
    147       CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    148       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    149      
     108    ENDDO
     109    !
     110    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
     111    !
     112    DO l = llm -1 , 2 , -1
     113       !
     114       DO ij = 1, ngrid
     115          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
     116          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
     117          beta (ij,l)  =   p(ij,l  ) / dellta   
     118       ENDDO
     119    ENDDO
    150120
    151       RETURN
    152       END
     121    !  ***********************************************************************
     122    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
     123    !
     124    DO   ij   = 1, ngrid
     125       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
     126            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
     127    ENDDO
     128    !
     129    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
     130    !
     131    DO l = 2, llm
     132       DO   ij   = 1, ngrid
     133          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
     134       ENDDO
     135    ENDDO
     136
     137    if (present(pkf)) then
     138       !    calcul de pkf
     139       pkf = pk
     140       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     141    end if
     142
     143  END SUBROUTINE exner_hyb
     144
     145end module exner_hyb_m
     146

  • trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu_m.F90

    r1300 r1302  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  F. Forget , Y. Wanherdrick
    7 c P.Le Van  , Fr. Hourdin  .
    8 c    ..........
    9 c
    10 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    11 c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    12 c
    13 c   ************************************************************************
    14 c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
    15 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    16 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    17 c   ************************************************************************
    18 c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    19 c    la pression et la fonction d'Exner  au  sol  .
    20 c
    21 c     WARNING : CECI est une version speciale de exner_hyb originale
    22 c               Utilise dans la version martienne pour pouvoir
    23 c               tourner avec des coordonnees verticales complexe
    24 c              => Il ne verifie PAS la condition la proportionalite en
    25 c              energie totale/ interne / potentielle (F.Forget 2001)
    26 c    ( voir note de Fr.Hourdin )  ,
    27 c
    28       IMPLICIT NONE
    29 c
    30 #include "dimensions.h"
    31 #include "paramet.h"
    32 #include "comconst.h"
    33 #include "comgeom.h"
    34 #include "comvert.h"
    35 #include "serre.h"
     1module exner_milieu_m
    362
    37       INTEGER  ngrid
    38       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    39       REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
     3  IMPLICIT NONE
    404
    41 c    .... variables locales   ...
     5contains
    426
    43       INTEGER l, ij
    44       REAL dum1
     7  SUBROUTINE  exner_milieu ( ngrid, ps, p, pks, pk, pkf )
     8    !
     9    !     Auteurs :  F. Forget , Y. Wanherdrick
     10    ! P.Le Van  , Fr. Hourdin  .
     11    !    ..........
     12    !
     13    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     14    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     15    !
     16    !   ************************************************************************
     17    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     18    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     19    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     20    !   ************************************************************************
     21    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     22    !    la pression et la fonction d'Exner  au  sol  .
     23    !
     24    !     WARNING : CECI est une version speciale de exner_hyb originale
     25    !               Utilise dans la version martienne pour pouvoir
     26    !               tourner avec des coordonnees verticales complexe
     27    !              => Il ne verifie PAS la condition la proportionalite en
     28    !              energie totale/ interne / potentielle (F.Forget 2001)
     29    !    ( voir note de Fr.Hourdin )  ,
     30    !
     31    !
     32    include "dimensions.h"
     33    include "paramet.h"
     34    include "comconst.h"
     35    include "comgeom.h"
     36    include "comvert.h"
     37    include "serre.h"
    4538
    46       REAL ppn(iim),pps(iim)
    47       REAL xpn, xps
    48       REAL SSUM
    49       EXTERNAL SSUM
    50       logical,save :: firstcall=.true.
    51       character(len=*),parameter :: modname="exner_milieu"
     39    INTEGER  ngrid
     40    REAL p(ngrid,llmp1),pk(ngrid,llm)
     41    real, optional:: pkf(ngrid,llm)
     42    REAL ps(ngrid),pks(ngrid)
    5243
    53       ! Sanity check
    54       if (firstcall) then
    55         ! sanity checks for Shallow Water case (1 vertical layer)
    56         if (llm.eq.1) then
     44    !    .... variables locales   ...
     45
     46    INTEGER l, ij
     47    REAL dum1
     48
     49    logical,save :: firstcall=.true.
     50    character(len=*),parameter :: modname="exner_milieu"
     51
     52    ! Sanity check
     53    if (firstcall) then
     54       ! sanity checks for Shallow Water case (1 vertical layer)
     55       if (llm.eq.1) then
    5756          if (kappa.ne.1) then
    58             call abort_gcm(modname,
    59      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     57             call abort_gcm(modname, &
     58                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6059          endif
    6160          if (cpp.ne.r) then
    62             call abort_gcm(modname,
    63      &      "cpp!=r , but running in Shallow Water mode!!",42)
     61             call abort_gcm(modname, &
     62                  "cpp!=r , but running in Shallow Water mode!!",42)
    6463          endif
    65         endif ! of if (llm.eq.1)
     64       endif ! of if (llm.eq.1)
    6665
    67         firstcall=.false.
    68       endif ! of if (firstcall)
     66       firstcall=.false.
     67    endif ! of if (firstcall)
    6968
    70 !!!! Specific behaviour for Shallow Water (1 vertical layer) case:
    71       if (llm.eq.1) then
    72      
    73         ! Compute pks(:),pk(:),pkf(:)
    74        
    75         DO   ij  = 1, ngrid
    76           pks(ij) = (cpp/preff) * ps(ij) 
     69    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     70    if (llm.eq.1) then
     71
     72       ! Compute pks(:),pk(:),pkf(:)
     73
     74       DO   ij  = 1, ngrid
     75          pks(ij) = (cpp/preff) * ps(ij)
    7776          pk(ij,1) = .5*pks(ij)
    78         ENDDO
    79        
    80         CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    81         CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    82        
    83         ! our work is done, exit routine
    84         return
     77       ENDDO
    8578
    86       endif ! of if (llm.eq.1)
     79       if (present(pkf)) then
     80          pkf = pk
     81          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     82       end if
    8783
    88 !!!! General case:
     84       ! our work is done, exit routine
     85       return
     86    endif ! of if (llm.eq.1)
    8987
    90 c     -------------
    91 c     Calcul de pks
    92 c     -------------
    93    
    94       DO   ij  = 1, ngrid
    95         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    96       ENDDO
     88    ! General case:
    9789
    98       DO  ij   = 1, iim
    99         ppn(ij) = aire(   ij   ) * pks(  ij     )
    100         pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    101       ENDDO
    102       xpn      = SSUM(iim,ppn,1) /apoln
    103       xps      = SSUM(iim,pps,1) /apols
     90    !     -------------
     91    !     Calcul de pks
     92    !     -------------
    10493
    105       DO ij   = 1, iip1
    106         pks(   ij     )  =  xpn
    107         pks( ij+ip1jm )  =  xps
    108       ENDDO
    109 c
    110 c
    111 c    .... Calcul de pk  pour la couche l
    112 c    --------------------------------------------
    113 c
    114       dum1 = cpp * (2*preff)**(-kappa)
    115       DO l = 1, llm-1
    116         DO   ij   = 1, ngrid
    117          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
    118         ENDDO
    119       ENDDO
     94    DO   ij  = 1, ngrid
     95       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     96    ENDDO
    12097
    121 c    .... Calcul de pk  pour la couche l = llm ..
    122 c    (on met la meme distance (en log pression)  entre Pk(llm)
    123 c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
     98    !    .... Calcul de pk  pour la couche l
     99    !    --------------------------------------------
     100    !
     101    dum1 = cpp * (2*preff)**(-kappa)
     102    DO l = 1, llm-1
     103       DO   ij   = 1, ngrid
     104          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
     105       ENDDO
     106    ENDDO
    124107
    125       DO   ij   = 1, ngrid
    126          pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    127       ENDDO
     108    !    .... Calcul de pk  pour la couche l = llm ..
     109    !    (on met la meme distance (en log pression)  entre Pk(llm)
     110    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
    128111
     112    DO   ij   = 1, ngrid
     113       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
     114    ENDDO
    129115
    130 c    calcul de pkf
    131 c    -------------
    132       CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    133       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    134      
    135 c    EST-CE UTILE ?? : calcul de beta
    136 c    --------------------------------
    137       DO l = 2, llm
    138         DO   ij   = 1, ngrid
    139           beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
    140         ENDDO
    141       ENDDO
     116    if (present(pkf)) then
     117       !    calcul de pkf
     118       pkf = pk
     119       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     120    end if
    142121
    143       RETURN
    144       END
     122  END SUBROUTINE exner_milieu
     123
     124end module exner_milieu_m
     125

  • trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90

    r1300 r1302  
    1414#endif
    1515  USE Write_Field
     16  use exner_hyb_m, only: exner_hyb
     17  use exner_milieu_m, only: exner_milieu
    1618
    1719  !   Author:    Frederic Hourdin      original: 15/01/93
    … …  
    5456  REAL pks(ip1jmp1)                      ! exner au  sol
    5557  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    56   REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    5758  REAL phi(ip1jmp1,llm)                  ! geopotentiel
    5859  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
    … …  
    223224        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    224225        if (pressure_exner) then
    225           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    226         else
    227           call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
     226          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
     227        else
     228          call exner_milieu(ip1jmp1,ps,p,pks,pk)
    228229        endif
    229230        CALL massdair(p,masse)

  • trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90

    r1300 r1302  
    7373  if (disvert_type==1) then
    7474     ! standard case for Earth (automatic generation of levels)
    75      call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig, scaleheight)
     75     call disvert()
    7676  else if (disvert_type==2) then
    7777     ! standard case for planets (levels generated using z2sig.def file)

  • trunk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F

    r1300 r1302  
    22! $Header$
    33!
    4 c
    5 c
     4!
     5!
    66
    77      subroutine q_sat(np,temp,pres,qsat)
    8 c
     8!
    99      IMPLICIT none
    10 c======================================================================
    11 c Autheur(s): Z.X. Li (LMD/CNRS)
    12 c  reecriture vectorisee par F. Hourdin.
    13 c Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
    14 c======================================================================
    15 c Arguments:
    16 c kelvin---input-R: temperature en Kelvin
    17 c millibar--input-R: pression en mb
    18 c
    19 c q_sat----output-R: vapeur d'eau saturante en kg/kg
    20 c======================================================================
    21 c
     10!======================================================================
     11! Autheur(s): Z.X. Li (LMD/CNRS)
     12!  reecriture vectorisee par F. Hourdin.
     13! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
     14!======================================================================
     15! Arguments:
     16! kelvin---input-R: temperature en Kelvin
     17! millibar--input-R: pression en mb
     18!
     19! q_sat----output-R: vapeur d'eau saturante en kg/kg
     20!======================================================================
     21!
    2222      integer np
    2323      REAL temp(np),pres(np),qsat(np)
    24 c
     24!
    2525      REAL r2es
    2626      PARAMETER (r2es=611.14 *18.0153/28.9644)
    27 c
     27!
    2828      REAL r3les, r3ies, r3es
    2929      PARAMETER (R3LES=17.269)
    3030      PARAMETER (R3IES=21.875)
    31 c
     31!
    3232      REAL r4les, r4ies, r4es
    3333      PARAMETER (R4LES=35.86)
    3434      PARAMETER (R4IES=7.66)
    35 c
     35!
    3636      REAL rtt
    3737      PARAMETER (rtt=273.16)
    38 c
     38!
    3939      REAL retv
    4040      PARAMETER (retv=28.9644/18.0153 - 1.0)
    … …  
    4242      real zqsat
    4343      integer ip
    44 c
    45 C     ------------------------------------------------------------------
    46 c
    47 c
     44!
     45!     ------------------------------------------------------------------
     46!
     47!
    4848
    4949      do ip=1,np
    5050
    51 c      write(*,*)'kelvin,millibar=',kelvin,millibar
    52 c       write(*,*)'temp,pres=',temp(ip),pres(ip)
    53 c
     51!      write(*,*)'kelvin,millibar=',kelvin,millibar
     52!       write(*,*)'temp,pres=',temp(ip),pres(ip)
     53!
    5454         IF (temp(ip) .LE. rtt) THEN
    5555            r3es = r3ies
    … …  
    5959            r4es = r4les
    6060         ENDIF
    61 c
     61!
    6262         zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
    6363         zqsat=MIN(0.5,ZQSAT)
    6464         zqsat=zqsat/(1.-retv *zqsat)
    65 c
     65!
    6666         qsat(ip)= zqsat
    67 c      write(*,*)'qsat=',qsat(ip)
     67!      write(*,*)'qsat=',qsat(ip)
    6868
    6969      enddo
    70 c
     70!
    7171      RETURN
    7272      END
     73
Note: See TracChangeset for help on using the changeset viewer.