Changeset 2868 for LMDZ5/trunk/libf


Ignore:
Timestamp:
May 2, 2017, 1:16:06 PM (8 years ago)
Author:
fhourdin
Message:

Option pour les fonctions de stabilité pour les coefficients
d'échanges en surface (Etienne Vignon)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cdrag.F90

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