Ignore:
Timestamp:
May 20, 2017, 9:41:16 AM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2865:2885 into testing branch

Location:
LMDZ5/branches/testing
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/cdrag.F90

    r2787 r2886  
    22!$Id$
    33!
    4  SUBROUTINE cdrag( knon,  nsrf,   &
     4 SUBROUTINE cdrag(knon,  nsrf,   &
    55     speed, t1,    q1,    zgeop1, &
    66     psol,  tsurf, qsurf, z0m, z0h,  &
    7      pcfm,  pcfh,  zri,   pref )
     7     cdm,  cdh,  zri,   pref)
    88
    99  USE dimphy
    1010  USE indice_sol_mod
    1111  USE print_control_mod, ONLY: lunout, prt_level
     12  USE ioipsl_getin_p_mod, ONLY : getin_p
     13
    1214  IMPLICIT NONE
    1315! ================================================================= c
    1416!
    1517! Objet : calcul des cdrags pour le moment (pcfm) et
    16 !         les flux de chaleur sensible et latente (pcfh).   
    17 !
    18 ! Modified histroy:
    19 !   27-Jan-2014: richardson number inconsistant between
    20 !   coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine
    21 !   by merging coefcdrag and clcdrag.
     18!         les flux de chaleur sensible et latente (pcfh) d'apr??s
     19!         Louis 1982, Louis 1979, King et al 2001
     20!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
     21!         et 1982 pour les cas instables
     22!
     23! Modified history:
     24!  writting on the 20/05/2016
     25!  modified on the 13/12/2016 to be adapted to LMDZ6
    2226!
    2327! References:
    … …  
    2832!     parametrization, November 1981, ECMWF, Reading, England.
    2933!     Page: 19. Equations in Table 1.
     34!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
     35!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
     36!    Boundary-Layer Meteorology 72: 331-344
    3037!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. 
    3138!     European Centre for Medium-Range Weather Forecasts.
    … …  
    3441!     model to the parameterization of evaporation from the tropical oceans. J.
    3542!     Climate, 5:418-434.
     43!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
     44!   to surface and boundary-layer flux parametrizations
     45!   QJRMS, 127, pp 779-794
    3646!
    3747! ================================================================= c
    38 !
    39 ! knon----input-I- nombre de points pour un type de surface
    40 ! nsrf----input-I- indice pour le type de surface; voir indicesol.h
    41 ! speed---input-R- module du vent au 1er niveau du modele
    42 ! t1------input-R- temperature de l'air au 1er niveau du modele
    43 ! q1------input-R- humidite de l'air au 1er niveau du modele
    44 ! zgeop---input-R- geopotentiel au 1er niveau du modele
    45 ! psol----input-R- pression au sol
    46 ! tsurf---input-R- temperature de l'air a la surface
    47 ! qsurf---input-R- humidite de l'air a la surface
    48 ! z0m, z0h---input-R- rugosite
    49 !! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015,
    50 !! u1------input-R- vent zonal au 1er niveau du modele
    51 !! v1------input-R- vent meridien au 1er niveau du modele
    52 !
    53 ! pcfm---output-R- cdrag pour le moment
    54 ! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
    55 ! zri----output-R- Richardson number
    56 ! pref---output-R- pression au niveau zgeop/RG
    57 !
    58 ! Parameters:
    59 ! ckap-----Karman constant
    60 ! cb,cc,cd-parameters in Louis et al., 1982
    6148! ================================================================= c
    62 !
    63 !
     49! On choisit le couple de fonctions de correction avec deux flags:
     50! Un pour les cas instables, un autre pour les cas stables
     51!
     52! iflag_corr_insta:
     53!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
     54!                   2: Louis 1982
     55!                   3: Laurent Li
     56!
     57! iflag_corr_sta:
     58!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
     59!                   2: Louis 1982
     60!                   3: Laurent Li
     61!                   4: King  2001 (SHARP)
     62!                   5: MO 1st order theory (allow collapse of turbulence)
     63!           
     64!
     65!*****************************************************************
    6466! Parametres d'entree
    6567!*****************************************************************
    66   INTEGER, INTENT(IN)                      :: knon, nsrf
     68 
     69  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
    6770  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
    6871  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
    69   REAL, DIMENSION(klon), INTENT(IN)        :: psol  ! pression au sol
    70   REAL, DIMENSION(klon), INTENT(IN)        :: t1    ! temperature at 1st level
    71   REAL, DIMENSION(klon), INTENT(IN)        :: q1    ! humidity at 1st level
    7272  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
    7373  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
    7474  REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m)
    75 !  paprs, pplay u1, v1: to be deleted
    76 !  they were in the old clcdrag. Fuxing WANG, 04/03/2015
    77 !  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
    78 !  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
    79 !  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1
     75  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K)
     76  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
     77  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
     78
     79
    8080
    8181! Parametres de sortie
    8282!******************************************************************
    83   REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm  ! Drag coefficient for heat flux
    84   REAL, DIMENSION(klon), INTENT(OUT)       :: pcfh  ! Drag coefficient for momentum
     83  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm  ! Drag coefficient for heat flux
     84  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh  ! Drag coefficient for momentum
    8585  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
    8686  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
    8787
    88 ! Parametres local
    89   INTEGER                                  :: ng_q1    ! Number of grids that q1 < 0.0
    90   INTEGER                                  :: ng_qsurf ! Number of grids that qsurf < 0.0
    91 !  zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015
    92 !  REAL, DIMENSION(klon)                    :: zgeop1! geopotentiel au 1er niveau du modele
    93 !  REAL, DIMENSION(klon)                    :: psol  ! pression au sol
    94 !
    95 ! ================================================================= c
    96 !
     88! Variables Locales
     89!******************************************************************
     90 
     91
    9792  INCLUDE "YOMCST.h"
    9893  INCLUDE "YOETHF.h"
    99 !  INCLUDE "indicesol.h"
    10094  INCLUDE "clesphys.h"
    101 !
    102 ! Quelques constantes et options:
    103 !!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
    104   REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2
    105 !
    106 ! Variables locales :
    107   INTEGER               :: i
    108   REAL                  :: zdu2, ztsolv
    109   REAL                  :: ztvd, zscf
    110   REAL                  :: zucf, zcr
    111   REAL                  :: friv, frih
    112   REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum
    113   REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux
    114   LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
    115   REAL, DIMENSION(klon) :: zcdn_m, zcdn_h         ! Drag coefficient in neutral conditions
     95
     96
     97  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
     98  REAL                   CEPDU2
     99  REAL                   ALPHA
     100  REAL                   CB,CC,CD,C2,C3
     101  REAL                   MU, CM, CH, B, CMstar, CHstar
     102  REAL                   PM, PH, BPRIME
     103  REAL                   C
     104  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
     105  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
     106  INTEGER                i
     107  REAL                   zdu2, ztsolv
     108  REAL                   ztvd, zscf
     109  REAL                   zucf, zcr
     110  REAL                   friv, frih
     111  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
     112  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
    116113  REAL zzzcd
    117 !
    118 ! Fonctions thermodynamiques et fonctions d'instabilite
    119   REAL                  :: fsta, fins, x
    120   fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
    121   fins(x) = SQRT(1.0-18.0*x)
     114
     115  LOGICAL, SAVE :: firstcall = .TRUE.
     116  !$OMP THREADPRIVATE(firstcall)
     117  INTEGER, SAVE :: iflag_corr_sta
     118  !$OMP THREADPRIVATE(iflag_corr_sta)
     119  INTEGER, SAVE :: iflag_corr_insta
     120  !$OMP THREADPRIVATE(iflag_corr_insta)
     121
     122!===================================================================c
     123! Valeurs numeriques des constantes
     124!===================================================================c
     125
     126
     127! Minimum du carre du vent
     128
     129 CEPDU2 = (0.1)**2
     130
     131! Louis 1982
     132
     133  CB=5.0
     134  CC=5.0
     135  CD=5.0
     136
     137
     138! King 2001
     139
     140  C2=0.25
     141  C3=0.0625
     142
     143 
     144! Louis 1979
     145
     146  BPRIME=4.7
     147  B=9.4
     148 
     149
     150!MO
     151
     152  ALPHA=5.0
     153 
    122154
    123155! ================================================================= c
    124 ! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1
    125 ! (le geopotentiel du premier couche de modele).
    126 ! zgeop1 is an input ivariable in this subroutine.
    127 !  DO i = 1, knon
    128 !     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
    129 !          * (paprs(i,1)-pplay(i,1))
    130 !  END DO
    131 ! ================================================================= c
    132 !
     156! Tests avant de commencer
    133157! Fuxing WANG, 04/03/2015
    134158! To check if there are negative q1, qsurf values.
     159!====================================================================c
    135160  ng_q1 = 0     ! Initialization
    136161  ng_qsurf = 0  ! Initialization
    … …  
    154179  ENDIF
    155180
    156 ! Calculer le frottement au sol (Cdrag)
    157   DO i = 1, knon
    158 !------------------------------------------------------------
    159 ! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015,
    160 !     zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2)
    161 !------------------------------------------------------------
     181
     182
     183!=============================================================================c
     184! Calcul du cdrag
     185!=============================================================================c
     186
     187! On choisit les fonctions de stabilite utilisees au premier appel
     188!**************************************************************************
     189  IF (firstcall) THEN
     190   iflag_corr_sta=2
     191   iflag_corr_insta=2
     192 
     193   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
     194   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
     195
     196   firstcall = .FALSE.
     197 ENDIF
     198
     199!xxxxxxxxxxxxxxxxxxxxxxx
     200  DO i = 1, knon        ! Boucle sur l'horizontal
     201!xxxxxxxxxxxxxxxxxxxxxxx
     202
     203
     204! calculs preliminaires:
     205!***********************
     206     
     207
    162208     zdu2 = MAX(CEPDU2, speed(i)**2)
    163 !     psol(i) = paprs(i,1)
    164209     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
    165                  (1.+ RETV * max(q1(i),0.0))))  ! negative q1 set to zero
    166 !------------ the old calculations in clcdrag----------------
    167 !     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i))
    168 !     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
    169 !          *(1.+RETV*q1(i))
    170 !------------------------------------------------------------
    171 ! Fuxing WANG, 04/03/2015, in this revised version,
    172 ! the negative qsurf and q1 are set to zero (as in coefcdrag)
    173      ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero
    174      ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
    175           *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero
     210                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
     211     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
     212     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
     213          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
    176214     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
    177215
    178216
    179 ! Coefficients CD neutres pour m et h : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h))
     217! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
     218!********************************************************************
     219
    180220     zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
    181      zcdn_m(i) = zzzcd*zzzcd
    182      zcdn_h(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
    183 
    184      IF (zri(i) .GT. 0.) THEN      ! situation stable
     221     cdmn(i) = zzzcd*zzzcd
     222     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
     223
     224
     225! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
     226!*******************************************************
     227
     228!''''''''''''''
     229! Cas instables
     230!''''''''''''''
     231
     232 IF (zri(i) .LT. 0.) THEN   
     233
     234
     235        SELECT CASE (iflag_corr_insta)
     236   
     237        CASE (1) ! Louis 1979 + Mascart 1995
     238
     239           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
     240           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
     241           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
     242           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
     243           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
     244           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     245            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
     246            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
     247           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     248            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
     249            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
     250         
     251
     252
     253
     254           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
     255           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
     256
     257        CASE (2) ! Louis 1982
     258
     259           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     260                *(1.0+zgeop1(i)/(RG*z0m(i)))))
     261           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     262           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     263
     264
     265        CASE (3) ! Laurent Li
     266
     267           
     268           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     269           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     270
     271
     272
     273         CASE default ! Louis 1982
     274
     275           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     276                *(1.0+zgeop1(i)/(RG*z0m(i)))))
     277           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     278           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     279
     280
     281         END SELECT
     282
     283
     284
     285! Calcul des drags
     286
     287
     288       cdm(i)=cdmn(i)*FM(i)
     289       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     290
     291
     292! Traitement particulier des cas oceaniques
     293! on applique Miller et al 1992 en l'absence de gustiness
     294
     295  IF (nsrf == is_oce) THEN
     296!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     297
     298        IF(iflag_gusts==0) THEN
     299           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
     300           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
     301        ENDIF
     302
     303
     304        cdm(i)=MIN(cdm(i),cdmmax)
     305        cdh(i)=MIN(cdh(i),cdhmax)
     306
     307  END IF
     308
     309
     310
     311 ELSE                         
     312
     313!'''''''''''''''
     314! Cas stables :
     315!'''''''''''''''
    185316        zri(i) = MIN(20.,zri(i))
    186         IF (.NOT.zxli) THEN
     317
     318       SELECT CASE (iflag_corr_sta)
     319   
     320        CASE (1) ! Louis 1979 + Mascart 1995
     321   
     322           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
     323           FH(i)=FM(i)
     324
     325
     326        CASE (2) ! Louis 1982
     327           
    187328           zscf = SQRT(1.+CD*ABS(zri(i)))
    188            friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
    189            zcfm1(i) = zcdn_m(i) * friv
    190            frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
    191 !!$ PB     zcfh1(i) = zcdn(i) * frih
    192 !!$ PB     zcfh1(i) = f_cdrag_stable * zcdn(i) * frih
    193            zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih
    194            IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih
    195 !!$ PB
    196            pcfm(i) = zcfm1(i)
    197            pcfh(i) = zcfh1(i)
    198         ELSE
    199            pcfm(i) = zcdn_m(i)* fsta(zri(i))
    200            pcfh(i) = zcdn_h(i)* fsta(zri(i))
    201         ENDIF
    202      ELSE                          ! situation instable
    203         IF (.NOT.zxli) THEN
    204            zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) &
    205                 *(1.0+zgeop1(i)/(RG*z0m(i)))))
    206            zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
    207 !!$ PB     zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
    208            zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
    209            pcfm(i) = zcfm2(i)
    210            pcfh(i) = zcfh2(i)
    211         ELSE
    212            pcfm(i) = zcdn_m(i)* fins(zri(i))
    213            pcfh(i) = zcdn_h(i)* fins(zri(i))
    214         ENDIF
    215         IF(iflag_gusts==0) THEN
    216 ! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active
    217            zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
    218            IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25)
    219         ENDIF
    220      ENDIF
    221   END DO
    222 
     329           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     330           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     331
     332
     333        CASE (3) ! Laurent Li
     334   
     335           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
     336           FH(i)=FM(i)
     337
     338
     339        CASE (4)  ! King 2001
     340         
     341           if (zri(i) .LT. C2/2.) then
     342           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
     343           FH(i)=  FM(i)
     344
     345
     346           else
     347           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
     348           FH(i)= FM(i)
     349           endif
     350
     351
     352        CASE (5) ! MO
     353   
     354          if (zri(i) .LT. 1./alpha) then
     355
     356           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
     357           FH(i)=FM(i)
     358
     359           else
     360
     361
     362           FM(i)=MAX(1E-7,f_ri_cd_min)
     363           FH(i)=FM(i)
     364
     365          endif
     366
     367
     368
     369
     370
     371        CASE default ! Louis 1982
     372
     373           zscf = SQRT(1.+CD*ABS(zri(i)))
     374           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     375           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     376
     377
     378
     379   END SELECT
     380
     381! Calcul des drags
     382
     383
     384       cdm(i)=cdmn(i)*FM(i)
     385       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     386
     387       IF(nsrf.EQ.is_oce) THEN
     388
     389        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     390        cdm(i)=MIN(cdm(i),cdmmax)
     391        cdh(i)=MIN(cdh(i),cdhmax)
     392
     393      ENDIF
     394
     395
     396 ENDIF
     397
     398
     399
     400
     401!xxxxxxxxxxx
     402  END DO   !  Fin de la boucle sur l'horizontal
     403!xxxxxxxxxxx
    223404! ================================================================= c
    224405     
    225   ! IM cf JLD : on seuille cdrag_m et cdrag_h
    226   IF (nsrf == is_oce) THEN
    227      DO i=1,knon
    228         pcfm(i)=MIN(pcfm(i),cdmmax)
    229         pcfh(i)=MIN(pcfh(i),cdhmax)
    230      END DO
    231   END IF
     406 
    232407
    233408END SUBROUTINE cdrag
     409
     410
     411!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  • LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90

    r2839 r2886  
    182182      CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
    183183   enddo
    184      CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
    185 #endif
    186 
     184#endif
     185
     186   CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
    187187   do k=1,PARASOL_NREFL
    188188    do ip=1, Npoints
    189      if (stlidar%cldlayer(ip,4).gt.0.01) then
     189     if (stlidar%cldlayer(ip,4).gt.0.01.and.stlidar%parasolrefl(ip,k).ne.missing_val) then
    190190       parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ &
    191191                            stlidar%cldlayer(ip,4)
    … …  
    196196     endif
    197197    enddo
     198   
    198199   enddo
    199200   CALL histwrite3d_cosp(o_Ncrefl,Ncref,nvertp)

  • LMDZ5/branches/testing/libf/phylmd/cpl_mod.F90

    r2546 r2886  
    5252  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp
    5353  !$OMP THREADPRIVATE(cpl_windsp)
     54  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_sens_rain, cpl_sens_snow
     55  !$OMP THREADPRIVATE(cpl_sens_rain, cpl_sens_snow)
    5456  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_taumod
    5557  !$OMP THREADPRIVATE(cpl_taumod)
    … …  
    9092  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp2D
    9193  !$OMP THREADPRIVATE(cpl_windsp2D)
     94  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: cpl_sens_rain2D, cpl_sens_snow2D
     95  !$OMP THREADPRIVATE(cpl_sens_rain2D, cpl_sens_snow2D)
    9296  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_atm_co22D
    9397  !$OMP THREADPRIVATE(cpl_atm_co22D)
    … …  
    169173    sum_error = sum_error + error
    170174    ALLOCATE(cpl_taumod(klon,2), stat = error)
     175    sum_error = sum_error + error
     176    ALLOCATE(cpl_sens_rain(klon,2), stat = error)
     177    sum_error = sum_error + error
     178    ALLOCATE(cpl_sens_snow(klon,2), stat = error)
    171179    sum_error = sum_error + error
    172180    ALLOCATE(cpl_rriv2D(nbp_lon,jj_nb), stat=error)
    … …  
    531539  SUBROUTINE cpl_send_ocean_fields(itime, knon, knindex, &
    532540       swdown, lwdown, fluxlat, fluxsens, &
    533        precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp)
     541       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp,&
     542       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
    534543!
    535544! This subroutine cumulates some fields for each time-step during a coupling
    … …  
    552561    REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder, albsol
    553562    REAL, DIMENSION(klon), INTENT(IN)       :: taux, tauy, windsp
     563    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
     564    REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
    554565
    555566! Local variables
    … …  
    583594       cpl_tauy(1:knon,cpl_index) = 0.0
    584595       cpl_windsp(1:knon,cpl_index) = 0.0
     596       cpl_sens_rain(1:knon,cpl_index) = 0.0
     597       cpl_sens_snow(1:knon,cpl_index) = 0.0
    585598       cpl_taumod(1:knon,cpl_index) = 0.0
    586599       IF (carbon_cycle_cpl) cpl_atm_co2(1:knon,cpl_index) = 0.0
    … …  
    614627       cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) + &
    615628            windsp(ig)      / REAL(nexca)
     629       cpl_sens_rain(ig,cpl_index) = cpl_sens_rain(ig,cpl_index) + &
     630            sens_prec_liq(ig)      / REAL(nexca)
     631       cpl_sens_snow(ig,cpl_index) = cpl_sens_snow(ig,cpl_index) + &
     632            sens_prec_sol(ig)      / REAL(nexca)
    616633       cpl_taumod(ig,cpl_index) =   cpl_taumod(ig,cpl_index) + &
    617634          SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL (nexca)
    … …  
    654671          sum_error = sum_error + error
    655672          ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
     673          sum_error = sum_error + error
     674          ALLOCATE(cpl_sens_rain2D(nbp_lon,jj_nb,2), stat=error)
     675          sum_error = sum_error + error
     676          ALLOCATE(cpl_sens_snow2D(nbp_lon,jj_nb,2), stat=error)
    656677          sum_error = sum_error + error
    657678          ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
    … …  
    706727            knon, knindex)
    707728
     729       CALL gath2cpl(cpl_sens_rain(:,cpl_index), cpl_sens_rain2D(:,:,cpl_index), &
     730            knon, knindex)
     731
     732       CALL gath2cpl(cpl_sens_snow(:,cpl_index), cpl_sens_snow2D(:,:,cpl_index), &
     733            knon, knindex)
     734
    708735       CALL gath2cpl(cpl_taumod(:,cpl_index), cpl_taumod2D(:,:,cpl_index), &
    709736            knon, knindex)
    … …  
    722749       pctsrf, lafin, rlon, rlat, &
    723750       swdown, lwdown, fluxlat, fluxsens, &
    724        precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy)
     751       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy,&
     752       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
    725753!
    726754! This subroutine cumulates some fields for each time-step during a coupling
    … …  
    746774    REAL, DIMENSION(klon), INTENT(IN)       :: albsol, taux, tauy
    747775    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
     776    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
     777    REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
    748778    LOGICAL, INTENT(IN)                     :: lafin
    749779
    … …  
    778808       cpl_taux(1:knon,cpl_index) = 0.0
    779809       cpl_tauy(1:knon,cpl_index) = 0.0
     810       cpl_sens_rain(1:knon,cpl_index) = 0.0
     811       cpl_sens_snow(1:knon,cpl_index) = 0.0
    780812       cpl_taumod(1:knon,cpl_index) = 0.0
    781813    ENDIF
    … …  
    806838       cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
    807839            tauy(ig)        / REAL(nexca)     
     840       cpl_sens_rain(ig,cpl_index) = cpl_sens_rain(ig,cpl_index) + &
     841            sens_prec_liq(ig)      / REAL(nexca)
     842       cpl_sens_snow(ig,cpl_index) = cpl_sens_snow(ig,cpl_index) + &
     843            sens_prec_sol(ig)      / REAL(nexca)
    808844       cpl_taumod(ig,cpl_index) = cpl_taumod(ig,cpl_index) + &
    809845            SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL(nexca)
    … …  
    839875          sum_error = sum_error + error
    840876          ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
     877          sum_error = sum_error + error
     878          ALLOCATE(cpl_sens_rain2D(nbp_lon,jj_nb,2), stat=error)
     879          sum_error = sum_error + error
     880          ALLOCATE(cpl_sens_snow2D(nbp_lon,jj_nb,2), stat=error)
    841881          sum_error = sum_error + error
    842882          ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
    … …  
    889929
    890930       CALL gath2cpl(cpl_tauy(:,cpl_index), cpl_tauy2D(:,:,cpl_index), &
     931            knon, knindex)
     932
     933       CALL gath2cpl(cpl_sens_rain(:,cpl_index), cpl_sens_rain2D(:,:,cpl_index), &
     934            knon, knindex)
     935
     936       CALL gath2cpl(cpl_sens_snow(:,cpl_index), cpl_sens_snow2D(:,:,cpl_index), &
    891937            knon, knindex)
    892938
    … …  
    10781124    tab_flds(:,:,ids_nsfice) = cpl_nsol2D(:,:,2)
    10791125    tab_flds(:,:,ids_dflxdt) = cpl_fder2D(:,:,2)
     1126    tab_flds(:,:,ids_qraioc) = cpl_sens_rain2D(:,:,1)
     1127    tab_flds(:,:,ids_qsnooc) = cpl_sens_snow2D(:,:,1)
     1128    tab_flds(:,:,ids_qraiic) = cpl_sens_rain2D(:,:,2)
     1129    tab_flds(:,:,ids_qsnoic) = cpl_sens_snow2D(:,:,2)
    10801130   
    10811131    IF (version_ocean=='nemo') THEN
    … …  
    12791329    DEALLOCATE(cpl_taux2D, cpl_tauy2D, cpl_windsp2D, cpl_taumod2D, stat=error )
    12801330    sum_error = sum_error + error
     1331    DEALLOCATE(cpl_sens_rain2D, cpl_sens_snow2D, stat=error)
     1332    sum_error = sum_error + error
     1333
    12811334   
    12821335    IF (carbon_cycle_cpl) THEN

  • LMDZ5/branches/testing/libf/phylmd/ener_conserv.F90

    r2870 r2886  
    2424USE phys_local_var_mod, ONLY : d_t_eva,d_t_lsc,d_q_eva,d_q_lsc
    2525USE phys_output_var_mod, ONLY : bils_ec,bils_ech,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss
     26USE add_phys_tend_mod, ONLY : fl_cor_ebil
    2627
    2728IMPLICIT none
    … …  
    5960   DO k = 1, klev
    6061   DO i = 1, klon
    61       ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
    62       d_t_ec(i,k)=0.5/ZRCPD &
    63  &      *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
    64       ENDDO
    65       ENDDO
     62     IF (fl_cor_ebil .GT. 0) then
     63       ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
     64     ELSE
     65       ZRCPD = RCPD*(1.0+RVTMP2*pqn(i,k))
     66     ENDIF
     67     d_t_ec(i,k)=0.5/ZRCPD &
     68 &     *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
     69   ENDDO
     70   ENDDO
    6671!-jld ec_conser
    6772

  • LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90

    r2839 r2886  
    368368                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
    369369         end if
     370        ENDDO
     371     ELSE  ! IF(k.LE.klevm1)
     372        DO i = 1, klon
     373           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
     374           zmqc(i) = 0.
    370375        ENDDO
    371376     ENDIF ! end IF(k.LE.klevm1)

  • LMDZ5/branches/testing/libf/phylmd/oasis.F90

    r2546 r2886  
    5555  INTEGER, PARAMETER :: ids_atmco2 = 24
    5656  INTEGER, PARAMETER :: ids_taumod = 25
    57   INTEGER, PARAMETER :: maxsend    = 25  ! Maximum number of fields to send
     57  INTEGER, PARAMETER :: ids_qraioc = 26
     58  INTEGER, PARAMETER :: ids_qsnooc = 27
     59  INTEGER, PARAMETER :: ids_qraiic = 28
     60  INTEGER, PARAMETER :: ids_qsnoic = 29
     61  INTEGER, PARAMETER :: maxsend    = 29  ! Maximum number of fields to send
    5862 
    5963  ! Id for fields received from ocean
    … …  
    177181            infosend(ids_atmco2)%action = .TRUE. ; infosend(ids_atmco2)%name = 'COATMCO2'
    178182        ENDIF
     183        infosend(ids_qraioc)%action = .TRUE. ; infosend(ids_qraioc)%name = 'COQRAIOC'
     184        infosend(ids_qsnooc)%action = .TRUE. ; infosend(ids_qsnooc)%name = 'COQSNOOC'
     185        infosend(ids_qraiic)%action = .TRUE. ; infosend(ids_qraiic)%name = 'COQRAIIC'
     186        infosend(ids_qsnoic)%action = .TRUE. ; infosend(ids_qsnoic)%name = 'COQSNOIC'
    179187       
    180188    ELSE IF (version_ocean=='opa8') THEN

  • LMDZ5/branches/testing/libf/phylmd/ocean_cpl_mod.F90

    r2546 r2886  
    191191    CALL cpl_send_ocean_fields(itime, knon, knindex, &
    192192         swnet, lwnet, fluxlat, fluxsens, &
    193          precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1, windsp&
    194          )
    195 !         ,sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
     193         precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1, windsp,&
     194         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
    196195   
    197196
    … …  
    353352       pctsrf, lafin, rlon, rlat, &
    354353       swnet, lwnet, fluxlat, fluxsens, &
    355        precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1 &
    356        )
    357 !      ,sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
     354       precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1,&
     355       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
    358356
    359357 

  • LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90

    r2870 r2886  
    23262326!!!
    23272327       IF (iflag_split .eq.0) THEN
     2328        wake_dltke(:,:,nsrf) = 0.
    23282329        DO k = 1, klev
    23292330           DO j = 1, knon

  • LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90

    r2870 r2886  
    352352!!!
    353353!!!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     354      LOGICAL, SAVE, ALLOCATABLE :: ptconv(:,:)
     355      !$OMP THREADPRIVATE(ptconv)
    354356!>jyg+nrlmd
    355357  !
    … …  
    676678      ALLOCATE(cdragm_x(klon), cdragm_w(klon))
    677679      ALLOCATE(kh(klon), kh_x(klon), kh_w(klon))
     680!
     681      ALLOCATE(ptconv(klon,klev))
    678682!
    679683      ALLOCATE(wbeff(klon), convoccur(klon), zmax_th(klon))
    … …  
    934938      DEALLOCATE(kh, kh_x, kh_w)
    935939!
     940      DEALLOCATE(ptconv)
     941!
    936942      DEALLOCATE(wbeff, convoccur, zmax_th)
    937943      DEALLOCATE(zq2m, zt2m, weak_inversion)

  • LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90

    r2870 r2886  
    867867      (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
    868868
    869   TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_qSTDlevs     = (/                             &
    870       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q850', "Specific humidity 850hPa", &
    871       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    872       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q700', "Specific humidity 700hPa", &
    873       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    874       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q500', "Specific humidity 500hPa", &
    875       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    876       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q200', "Specific humidity 200hPa", &
    877       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    878       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q100', "Specific humidity 100hPa", &
    879       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    880       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q50', "Specific humidity 50hPa",  &
    881       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
    882       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q10', "Specific humidity 10hPa", &
    883       "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
     869  TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_qSTDlevs     = (/ &
     870       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q850', &
     871       "Specific humidity 850hPa", "kg/kg", &
     872       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     873       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q700', &
     874       "Specific humidity 700hPa", "kg/kg", &
     875       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     876       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q500', &
     877       "Specific humidity 500hPa", "kg/kg", &
     878       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     879       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q200', &
     880       "Specific humidity 200hPa", "kg/kg", &
     881       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     882       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q100', &
     883       "Specific humidity 100hPa", "kg/kg", &
     884       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     885       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q50', &
     886       "Specific humidity 50hPa", "kg/kg", &
     887       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
     888       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q10', &
     889       "Specific humidity 10hPa", "kg/kg", &
     890       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
    884891
    885892  TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_zSTDlevs   = (/                           &

  • LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90

    r2870 r2886  
    156156!!!       d_s_the, d_dens_the, &            ! due to thermals
    157157       !                                 
     158       ptconv, &
    158159       wbeff, convoccur, zmax_th, &
    159160       sens, flwp, fiwp,  &
    … …  
    599600    LOGICAL,SAVE :: ok_adjwk=.FALSE.
    600601    !$OMP THREADPRIVATE(ok_adjwk)
     602    INTEGER,SAVE :: iflag_adjwk=0            !jyg
     603    !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
    601604    REAL,SAVE :: oliqmax=999.,oicemax=999.
    602605    !$OMP THREADPRIVATE(oliqmax,oicemax)
    … …  
    876879    save iflag_cld_th
    877880    !$OMP THREADPRIVATE(iflag_cld_th)
    878     logical ptconv(klon,klev)
     881!IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
    879882    !IM cf. AM 081204 BEG
    880883    logical ptconvth(klon,klev)
    … …  
    12011204       CALL suphel ! initialiser constantes et parametres phys.
    12021205       CALL getin_p('random_notrig_max',random_notrig_max)
    1203        CALL getin_p('ok_adjwk',ok_adjwk)
     1206       CALL getin_p('ok_adjwk',ok_adjwk)
     1207       IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
     1208       ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
     1209                      ! 1 => convective adjustment but state variables are unchanged
     1210                      ! 2 => convective adjustment and state variables are changed
     1211       CALL getin_p('iflag_adjwk',iflag_adjwk)
    12041212       CALL getin_p('oliqmax',oliqmax)
    12051213       CALL getin_p('oicemax',oicemax)
    … …  
    24232431       ! after the call to the convective scheme.
    24242432       IF (iflag_wake>=1) then
    2425           IF (ok_adjwk) THEN
     2433          IF (iflag_adjwk >= 1) THEN
    24262434             limbas(:) = 1
    24272435             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
    … …  
    24412449                ENDDO
    24422450             ENDDO
    2443              CALL add_wake_tend &
     2451             IF (iflag_adjwk == 2) THEN
     2452               CALL add_wake_tend &
    24442453                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
    2445           ENDIF  ! (ok_adjwk)
     2454             ENDIF  ! (iflag_adjwk == 2)
     2455          ENDIF  ! (iflag_adjwk >= 1)
    24462456       ENDIF ! (iflag_wake>=1)
    24472457       !>jyg
    … …  
    25252535          !    Add the tendency due to the dry adjustment of the wake profile
    25262536          IF (iflag_wake>=1) THEN
    2527              DO k=1,klev
    2528                 DO i=1,klon
    2529                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
    2530                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
    2531                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
    2532                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
    2533                 ENDDO
    2534              ENDDO
    2535           ENDIF
     2537            IF (iflag_adjwk == 2) THEN
     2538              DO k=1,klev
     2539                 DO i=1,klon
     2540                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
     2541                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
     2542                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
     2543                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
     2544                 ENDDO
     2545              ENDDO
     2546            ENDIF  ! (iflag_adjwk = 2)
     2547          ENDIF   ! (iflag_wake>=1)
    25362548          !>jyg
    25372549          !

  • LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_init_140gp.F90

    r2870 r2886  
    8181!  Calculate lookup tables for functions needed in routine TAUMOL (TAUGB2)
    8282
     83! FH 2017/05/03
     84! Ce facteur de correction CORR2 est vraiment bizare parce qu'on
     85! impose 1. aux bornes,  en I=1 et I=200 mais la fonction
     86! CORE=( 1 - sqrt(i/im) ) / ( 1 - i/im ) = 1/ ( 1 + sqrt(i/im))
     87! vaut 1 en i=1 et 1/2 en i=im ...
     88
    8389CORR1(0) = 1.0_JPRB
    84 CORR1(400) = 1.0_JPRB
     90CORR1(200) = 1.0_JPRB
    8591CORR2(0) = 1.0_JPRB
    86 CORR2(400) = 1.0_JPRB
    87 DO I = 1,399
    88   Z_FP = 0.0025_JPRB*REAL(I)
     92CORR2(200) = 1.0_JPRB
     93DO I = 1,199
     94  Z_FP = 0.005_JPRB*REAL(I)
    8995  Z_RTFP = SQRT(Z_FP)
    9096  CORR1(I) = Z_RTFP/Z_FP

  • LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_taumol2.F90

    r1999 r2886  
    111111  IFP=MAX(0,IFP)
    112112
     113! FH 2017/05/02
     114! Modification parce qu'on avait un plantage sur un cas 1D.
     115! C'est evidemment une correction suspecte
     116  IF (IFP>200) THEN
     117      PRINT*,'WARNING IFP=',IFP,' 2.E2_JPRB*Z_FP+0.5_JPRB avec Z_FP=',Z_FP
     118      IFP=200
     119  ENDIF
     120
    113121  Z_FC00(I_LAY) = P_FAC00(I_LAY) * CORR2(IFP)
    114122  Z_FC10(I_LAY) = P_FAC10(I_LAY) * CORR2(IFP)
    … …  
    147155!---MI 981104       
    148156  IF (IFP <= 0) IFP=0
     157  IF (IFP>200) THEN
     158      PRINT*,'WARNING IFP=',IFP,' 2.E2_JPRB*Z_FP+0.5_JPRB avec Z_FP=',Z_FP
     159      IFP=200
     160  ENDIF
    149161
    150162  Z_FC00(I_LAY) = P_FAC00(I_LAY) * CORR2(IFP)

  • LMDZ5/branches/testing/libf/phylmd/rrtm/yoerrtbg2.F90

    r2870 r2886  
    1313!    -------------------------------------------------------------------
    1414
    15 REAL(KIND=JPRB) :: CORR1(0:400)
    16 REAL(KIND=JPRB) :: CORR2(0:400)
     15REAL(KIND=JPRB) :: CORR1(0:200)
     16REAL(KIND=JPRB) :: CORR2(0:200)
    1717
    1818!     -----------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.