Ignore:
Timestamp:
Aug 31, 2023, 2:58:34 PM (17 months ago)
Author:
Laurent Fairhead
Message:

Inclusion and merge of Olivier's modifications

Location:
LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd
Files:
9 added
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/cdrag_mod.F90

    r4481 r4661  
    1515 SUBROUTINE cdrag(knon,  nsrf,   &
    1616     speed, t1,    q1,    zgeop1, &
    17      psol, tsurf, qsurf, z0m, z0h,  &
     17     psol, pblh, tsurf, qsurf, z0m, z0h,  &
    1818     ri_in, iri_in, &
    19      cdm,  cdh,  zri,   pref)
     19     cdm,  cdh,  zri,   pref, prain, tsol , pat1)
    2020
    2121  USE dimphy
    … …  
    8181 
    8282  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
    83   REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
    84   REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
    85   REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
    86   REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
    87   REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m)
     83  REAL, DIMENSION(klon), INTENT(IN)        :: speed    ! module du vent au 1er niveau du modele
     84  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1   ! geopotentiel au 1er niveau du modele
     85  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf    ! Surface temperature (K)
     86  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf    ! Surface humidity (Kg/Kg)
     87  REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m, z0h ! Rugosity at surface (m)
    8888  REAL, DIMENSION(klon), INTENT(IN)        :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var.
    8989  INTEGER, INTENT(IN)                      :: iri_in! iflag to activate cdrag calculation using ri1
    90   REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K)
    91   REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
    92   REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
     90  REAL, DIMENSION(klon), INTENT(IN)        :: t1       ! Temperature au premier niveau (K)
     91  REAL, DIMENSION(klon), INTENT(IN)        :: q1       ! humidite specifique au premier niveau (kg/kg)
     92  REAL, DIMENSION(klon), INTENT(IN)        :: psol     ! pression au sol
     93
     94!------------------ Rajout pour COARE (OT2018) --------------------
     95  REAL, DIMENSION(klon), INTENT(INOUT)     :: pblh  !hauteur de CL
     96  REAL, DIMENSION(klon), INTENT(IN)        :: prain !rapport au precipitation
     97  REAL, DIMENSION(klon), INTENT(IN)        :: tsol  !SST imposé sur la surface oceanique
     98  REAL, DIMENSION(klon), INTENT(IN)        :: pat1  !pression premier lev
    9399
    94100
    … …  
    99105  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh   ! Drag coefficient for heat flux
    100106  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
    101   REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
     107  REAL, DIMENSION(klon), INTENT(INOUT)       :: pref  ! Pression au niveau zgeop/RG
    102108
    103109! Variables Locales
    … …  
    119125  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
    120126  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
    121   INTEGER                i
     127  INTEGER                i, k
    122128  REAL                   zdu2, ztsolv
    123129  REAL                   ztvd, zscf
    … …  
    129135  REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
    130136  REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl
     137
     138!------------------ Rajout (OT2018) --------------------
     139!------------------ Rajout pour les appelles BULK (OT) --------------------
     140  REAL, DIMENSION(klon,2) :: rugos_itm
     141  REAL, DIMENSION(klon,2) :: rugos_ith
     142  REAL, PARAMETER         :: tol_it_rugos=1.e-4
     143  REAL, DIMENSION(3)      :: coeffs
     144  LOGICAL                 :: mixte
     145  REAL                    :: z_0m
     146  REAL                    :: z_0h
     147  REAL, DIMENSION(klon)   :: cdmm
     148  REAL, DIMENSION(klon)   :: cdhh
     149
     150!--------- rajout pour verification --------------------
     151
     152  REAL, DIMENSION(klon) :: zcdn_h2
     153
     154
     155!------------------RAJOUT POUR ECUME -------------------
     156
     157  REAL, DIMENSION(klon) :: PQSAT
     158  REAL, DIMENSION(klon) :: PSFTH
     159  REAL, DIMENSION(klon) :: PFSTQ
     160  REAL, DIMENSION(klon) :: PUSTAR
     161  REAL, DIMENSION(klon) :: PCD      ! Drag coefficient for momentum
     162  REAL, DIMENSION(klon) :: PCDN     ! Drag coefficient for momentum
     163  REAL, DIMENSION(klon) :: PCH      ! Drag coefficient for momentum
     164  REAL, DIMENSION(klon) :: PCE      ! Drag coefficient for momentum
     165  REAL, DIMENSION(klon) :: PRI
     166  REAL, DIMENSION(klon) :: PRESA
     167  REAL, DIMENSION(klon) :: PSSS
     168
     169  LOGICAL               :: OPRECIP
     170  LOGICAL               :: OPWEBB
     171  LOGICAL               :: OPERTFLUX
     172  LOGICAL               :: LPRECIP
     173  LOGICAL               :: LPWG
     174
     175
    131176
    132177  LOGICAL, SAVE :: firstcall = .TRUE.
    … …  
    219264 ENDIF
    220265
     266!------------------ Rajout (OT2018) --------------------
     267!---------    Rajout pour itération sur rugosité     ----------------
     268  rugos_itm(:,2) = z0m
     269  rugos_itm(:,1) = 3.*tol_it_rugos*z0m
     270
     271  rugos_ith(:,2) = z0h !cp nouveau rugos_it
     272  rugos_ith(:,1) = 3.*tol_it_rugos*z0h
     273!--------------------------------------------------------------------
     274
    221275!xxxxxxxxxxxxxxxxxxxxxxx
    222276  DO i = 1, knon        ! Boucle sur l'horizontal
    … …  
    227281!***********************
    228282     
    229 
    230283     zdu2 = MAX(CEPDU2, speed(i)**2)
    231284     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
    232                  (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
     285          (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
    233286     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
    234      ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
     287     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
    235288          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
    236      zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     289     
     290!------------------ Rajout (OT2018) --------------------
     291!--------------   ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN    ------------------------     
    237292     IF (iri_in.EQ.1) THEN
    238293       zri(i) = ri_in(i)
    239294     ENDIF
    240295
     296     IF (nsrf == is_oce) THEN
     297       
     298!------------------  Pour Core 2 choix Core Pur ou Core Mixte  --------------------
     299       IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN
     300         IF(choix_bulk .eq. 2) THEN
     301           mixte = .false.
     302         ELSE
     303            mixte = .true.
     304         ENDIF
     305         call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),&
     306             zgeop1(i)/RG, zgeop1(i)/RG,  zgeop1(i)/RG,&
     307             psol(i),nit_bulk,mixte,&
     308             coeffs,z_0m,z_0h)
     309         cdmm(i) = coeffs(1)
     310         cdhh(i) = coeffs(2)
     311         cdm(i)=cdmm(i)
     312         cdh(i)=cdhh(i)
     313         write(*,*) "clc_core cd ch",cdmm(i),cdhh(i)
     314
     315!------------------                 Pour ECUME                 --------------------     
     316       ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN
     317         OPRECIP = .false.
     318         OPWEBB  = .false.
     319         OPERTFLUX = .false.
     320         IF (nsrf .eq. is_oce) THEN
     321           PSSS = 0.0
     322         ENDIF
     323         call ini_csts
     324         call ecumev6_flux( z_0m,t1(i),tsurf(i),&
     325             q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,&
     326             psol(i),pat1(i), OPRECIP, OPWEBB,&
     327             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,&
     328             PCH,PCE,PRI,PRESA,prain,&
     329             z_0h,OPERTFLUX,coeffs)
     330       
     331         cdmm(i) = coeffs(1)
     332         cdhh(i) = coeffs(2)
     333         cdm(i)=cdmm(i)
     334         cdh(i)=cdhh(i)
     335   
     336         write(*,*) "ecume cd ch",cdmm(i),cdhh(i)
     337       
     338!------------------                 Pour COARE CNRM                 --------------------
     339       ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN
     340         LPRECIP = .false.
     341         LPWG    = .false.
     342         call ini_csts
     343         call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i),  &
     344             sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, &
     345             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
     346             PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs)
     347         cdmm(i) = coeffs(1)
     348         cdhh(i) = coeffs(2)
     349         cdm(i)=cdmm(i)
     350         cdh(i)=cdhh(i)
     351         write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i)
     352
     353!------------------                 Pour COARE Maison                 --------------------
     354       ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN
     355         IF ( pblh(i) .eq. 0. ) THEN
     356           pblh(i) = 1500.
     357         ENDIF
     358         write(*,*) "debug size",size(coeffs)
     359         call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),&
     360                 t1(i),q1(i),&
     361                 zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,&
     362                 psol(i), pblh(i),&
     363                 nit_bulk,&
     364                 coeffs,z_0m,z_0h)
     365         cdmm(i) = coeffs(1)
     366         cdhh(i) = coeffs(3)
     367         cdm(i)=cdmm(i)
     368         cdh(i)=cdhh(i)
     369         write(*,*) "coare cd ch",cdmm(i),cdhh(i)
     370       ELSE
     371!------------------                 Pour La param LMDZ (ocean)              --------------------
     372         zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     373         IF (iri_in.EQ.1) THEN
     374           zri(i) = ri_in(i)
     375         ENDIF
     376       ENDIF
     377
     378!----------------------- Rajout des itérations --------------
     379       DO  k=1,nit_bulk
    241380
    242381! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
    243382!********************************************************************
    244 
    245      zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
    246      cdmn(i) = zzzcd*zzzcd
    247      cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
     383         zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2)))
     384         cdmn(i) = zzzcd*zzzcd
     385         cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2))))
     386
     387! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
     388!*******************************************************
     389         IF (zri(i) .LT. 0.) THEN   
     390           SELECT CASE (iflag_corr_insta)
     391             CASE (1) ! Louis 1979 + Mascart 1995
     392                  MU=LOG(MAX(z0m(i)/z0h(i),0.01))
     393                  CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
     394                  PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
     395                  CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
     396                  PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
     397                  CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     398                     & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
     399                     & * ((zgeop1(i)/(RG*z0h(i)))**PH)
     400                  CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     401                     & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
     402                     & * ((zgeop1(i)/(RG*z0m(i)))**PM)
     403                  FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
     404                  FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
     405             CASE (2) ! Louis 1982
     406                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     407                        *(1.0+zgeop1(i)/(RG*z0m(i)))))
     408                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     409                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     410             CASE (3) ! Laurent Li
     411                  FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     412                  FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     413             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     414                  sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
     415                  prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
     416                  FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
     417                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     418             CASE default ! Louis 1982
     419                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     420                         *(1.0+zgeop1(i)/(RG*z0m(i)))))
     421                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     422                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     423           END SELECT
     424! Calcul des drags
     425           cdmm(i)=cdmn(i)*FM(i)
     426           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     427! Traitement particulier des cas oceaniques
     428! on applique Miller et al 1992 en l'absence de gustiness
     429           IF (nsrf == is_oce) THEN
     430!            cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     431             IF (iflag_gusts==0) THEN
     432               zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
     433               cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
     434             ENDIF
     435             cdmm(i)=MIN(cdmm(i),cdmmax)
     436             cdhh(i)=MIN(cdhh(i),cdhmax)
     437             write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i)
     438           END IF
     439         ELSE                         
     440
     441!'''''''''''''''
     442! Cas stables :
     443!'''''''''''''''
     444           zri(i) = MIN(20.,zri(i))
     445           SELECT CASE (iflag_corr_sta)
     446             CASE (1) ! Louis 1979 + Mascart 1995
     447                  FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
     448                  FH(i)=FM(i)
     449             CASE (2) ! Louis 1982
     450                  zscf = SQRT(1.+CD*ABS(zri(i)))
     451                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     452                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     453             CASE (3) ! Laurent Li
     454                  FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
     455                  FH(i)=FM(i)
     456             CASE (4)  ! King 2001
     457                  IF (zri(i) .LT. C2/2.) THEN
     458                    FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
     459                    FH(i)=  FM(i)
     460                  ELSE
     461                    FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
     462                    FH(i)= FM(i)
     463                  ENDIF
     464             CASE (5) ! MO
     465                  if (zri(i) .LT. 1./alpha) then
     466                    FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
     467                    FH(i)=FM(i)
     468                  else
     469                    FM(i)=MAX(1E-7,f_ri_cd_min)
     470                    FH(i)=FM(i)
     471                  endif
     472             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     473                  sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
     474                  prandtl(i) = pr_neut + zri(i) * pr_slope
     475                  FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
     476                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     477             CASE default ! Louis 1982
     478                  zscf = SQRT(1.+CD*ABS(zri(i)))
     479                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     480                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     481           END SELECT
     482! Calcul des drags
     483           cdmm(i)=cdmn(i)*FM(i)
     484           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     485           IF (nsrf.EQ.is_oce) THEN
     486             cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     487             cdmm(i)=MIN(cdm(i),cdmmax)
     488             cdhh(i)=MIN(cdh(i),cdhmax)
     489           ENDIF
     490           rugos_itm(i,1) = rugos_itm(i,2)
     491           rugos_ith(i,1) = rugos_ith(i,2)
     492           rugos_itm(i,2) =  0.018*cdmm(i) * (speed(i))/RG  &
     493                              +  0.11*14e-6 / SQRT(cdmm(i) * zdu2)
     494
     495!---------- Version SEPARATION DES Z0 ----------------------
     496           IF (iflag_z0_oce==0) THEN
     497             rugos_ith(i,2) = rugos_itm(i,2)
     498           ELSE IF (iflag_z0_oce==1) THEN
     499             rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2)
     500           ENDIF
     501         ENDIF
     502         rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2))
     503         rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2))
     504       ENDDO
     505       IF (nsrf.EQ.is_oce) THEN
     506         cdm(i)=MIN(cdmm(i),cdmmax)
     507         cdh(i)=MIN(cdhh(i),cdhmax)
     508       ENDIF
     509       z0m = rugos_itm(:,2)
     510       z0h = rugos_ith(:,2)
     511     ELSE ! (nsrf == is_oce)
     512       zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     513       IF (iri_in.EQ.1) THEN
     514         zri(i) = ri_in(i)
     515       ENDIF
     516
     517! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
     518!********************************************************************
     519       zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
     520       cdmn(i) = zzzcd*zzzcd
     521       cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
    248522
    249523
    250524! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
    251525!*******************************************************
    252 
    253 
    254 
    255526!''''''''''''''
    256527! Cas instables
    257528!''''''''''''''
    258 
    259  IF (zri(i) .LT. 0.) THEN   
    260 
    261 
    262         SELECT CASE (iflag_corr_insta)
    263    
    264         CASE (1) ! Louis 1979 + Mascart 1995
    265 
    266            MU=LOG(MAX(z0m(i)/z0h(i),0.01))
    267            CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
    268            PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
    269            CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
    270            PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
    271            CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
    272             & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
    273             & * ((zgeop1(i)/(RG*z0h(i)))**PH)
    274            CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
    275             & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
    276             & * ((zgeop1(i)/(RG*z0m(i)))**PM)
    277          
    278 
    279 
    280 
    281            FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
    282            FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
    283 
    284         CASE (2) ! Louis 1982
    285 
    286            zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
    287                 *(1.0+zgeop1(i)/(RG*z0m(i)))))
    288            FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
    289            FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
    290 
    291 
    292         CASE (3) ! Laurent Li
    293 
    294            
    295            FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
    296            FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
    297 
    298          CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
    299          
    300            sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
    301            prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
    302            FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
    303            FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
    304 
    305          CASE default ! Louis 1982
    306 
    307            zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
    308                 *(1.0+zgeop1(i)/(RG*z0m(i)))))
    309            FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
    310            FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
    311 
    312 
     529       IF (zri(i) .LT. 0.) THEN   
     530         SELECT CASE (iflag_corr_insta)
     531           CASE (1) ! Louis 1979 + Mascart 1995
     532                MU=LOG(MAX(z0m(i)/z0h(i),0.01))
     533                CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
     534                PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
     535                CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
     536                PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
     537                CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     538                   & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
     539                   & * ((zgeop1(i)/(RG*z0h(i)))**PH)
     540                CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     541                   & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
     542                   & * ((zgeop1(i)/(RG*z0m(i)))**PM)
     543                FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
     544                FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
     545           CASE (2) ! Louis 1982
     546                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     547                       *(1.0+zgeop1(i)/(RG*z0m(i)))))
     548                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     549                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     550           CASE (3) ! Laurent Li
     551                FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     552                FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     553           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     554                 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
     555                 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
     556                 FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
     557                 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     558           CASE default ! Louis 1982
     559                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     560                      *(1.0+zgeop1(i)/(RG*z0m(i)))))
     561                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     562                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
    313563         END SELECT
    314 
    315 
    316 
    317564! Calcul des drags
    318 
    319 
    320        cdm(i)=cdmn(i)*FM(i)
    321        cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
    322 
    323 
    324 ! Traitement particulier des cas oceaniques
    325 ! on applique Miller et al 1992 en l'absence de gustiness
    326 
    327   IF (nsrf == is_oce) THEN
    328 !        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
    329 
    330         IF(iflag_gusts==0) THEN
    331            zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
    332            cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
    333         ENDIF
    334 
    335 
    336         cdm(i)=MIN(cdm(i),cdmmax)
    337         cdh(i)=MIN(cdh(i),cdhmax)
    338 
    339   END IF
    340 
    341 
    342 
    343  ELSE                         
    344 
     565         cdm(i)=cdmn(i)*FM(i)
     566         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     567       ELSE                         
    345568!'''''''''''''''
    346569! Cas stables :
    347570!'''''''''''''''
    348         zri(i) = MIN(20.,zri(i))
    349 
    350        SELECT CASE (iflag_corr_sta)
    351    
    352         CASE (1) ! Louis 1979 + Mascart 1995
    353    
    354            FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
    355            FH(i)=FM(i)
    356 
    357 
    358         CASE (2) ! Louis 1982
    359            
    360            zscf = SQRT(1.+CD*ABS(zri(i)))
    361            FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
    362            FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
    363 
    364 
    365         CASE (3) ! Laurent Li
    366    
    367            FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
    368            FH(i)=FM(i)
    369 
    370 
    371         CASE (4)  ! King 2001
    372          
    373            if (zri(i) .LT. C2/2.) then
    374            FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
    375            FH(i)=  FM(i)
    376 
    377 
    378            else
    379            FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
    380            FH(i)= FM(i)
    381            endif
    382 
    383 
    384         CASE (5) ! MO
    385    
    386           if (zri(i) .LT. 1./alpha) then
    387 
    388            FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
    389            FH(i)=FM(i)
    390 
    391            else
    392 
    393 
    394            FM(i)=MAX(1E-7,f_ri_cd_min)
    395            FH(i)=FM(i)
    396 
    397           endif
    398 
    399         CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
    400           sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
    401           prandtl(i) = pr_neut + zri(i) * pr_slope
    402           FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
    403           FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
    404 
    405         CASE default ! Louis 1982
    406 
    407            zscf = SQRT(1.+CD*ABS(zri(i)))
    408            FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
    409            FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
    410 
    411 
    412 
    413    END SELECT
    414 
     571         zri(i) = MIN(20.,zri(i))
     572         SELECT CASE (iflag_corr_sta)
     573           CASE (1) ! Louis 1979 + Mascart 1995
     574                FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
     575                FH(i)=FM(i)
     576           CASE (2) ! Louis 1982
     577                zscf = SQRT(1.+CD*ABS(zri(i)))
     578                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     579                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     580           CASE (3) ! Laurent Li
     581                FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
     582                FH(i)=FM(i)
     583           CASE (4)  ! King 2001
     584                if (zri(i) .LT. C2/2.) then
     585                  FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
     586                  FH(i)=  FM(i)
     587                else
     588                  FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
     589                  FH(i)= FM(i)
     590                endif
     591           CASE (5) ! MO
     592                if (zri(i) .LT. 1./alpha) then
     593                  FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
     594                  FH(i)=FM(i)
     595                else
     596                  FM(i)=MAX(1E-7,f_ri_cd_min)
     597                  FH(i)=FM(i)
     598                endif
     599           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     600                sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
     601                prandtl(i) = pr_neut + zri(i) * pr_slope
     602                FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
     603                FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     604           CASE default ! Louis 1982
     605                zscf = SQRT(1.+CD*ABS(zri(i)))
     606                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     607                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     608         END SELECT
    415609! Calcul des drags
    416 
    417 
    418        cdm(i)=cdmn(i)*FM(i)
    419        cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
    420 
    421        IF(nsrf.EQ.is_oce) THEN
    422 
    423         cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
    424         cdm(i)=MIN(cdm(i),cdmmax)
    425         cdh(i)=MIN(cdh(i),cdhmax)
    426 
    427       ENDIF
    428 
    429 
    430  ENDIF
    431 
    432 !xxxxxxxxxxx
     610         cdm(i)=cdmn(i)*FM(i)
     611         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     612       ENDIF
     613     ENDIF ! fin du if (nsrf == is_oce)
    433614  END DO   !  Fin de la boucle sur l'horizontal
    434 !xxxxxxxxxxx
    435 ! ================================================================= c
    436615     
    437616END SUBROUTINE cdrag

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/clesphys.h

    r4458 r4661  
    4242!IM seuils cdrm, cdrh
    4343       REAL cdmmax, cdhmax
     44!IM pour les params différentes Olivier Torres
     45       INTEGER choix_bulk, nit_bulk, kz0
    4446!IM param. stabilite s/ terres et en dehors
    4547       REAL ksta, ksta_ter, f_ri_cd_min
    … …  
    105107       COMMON/clesphys/                                                 &
    106108! REAL FIRST
     109! rajout choix_bulk et nit_bulk kz0 par Olivier Torres
    107110     &       co2_ppm, solaire                                           &
    108111     &     , RCO2, RCH4, RN2O, RCFC11, RCFC12                           &
    … …  
    132135     &     , ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad                &
    133136     &     , iflag_con, nbapp_cv, nbapp_wk                              &
     137     &     , choix_bulk, nit_bulk, kz0                                  &
    134138     &     , iflag_ener_conserv                                         &
    135139     &     , ok_suntime_rrtm                                            &

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/conf_phys_m.F90

    r4466 r4661  
    255255    LOGICAL, SAVE :: ok_new_lscp_omp
    256256    LOGICAL, SAVE :: ok_icefra_lscp_omp
     257    !rajout de choix_bulk et nit_bulk par Olivier Torres
     258    INTEGER,SAVE  :: choix_bulk_omp
     259    INTEGER,SAVE  :: nit_bulk_omp
     260    INTEGER,SAVE  :: kz0_omp
    257261
    258262
    … …  
    947951    nbapp_rad_omp = 12
    948952    CALL getin('nbapp_rad',nbapp_rad_omp)
     953
     954    !rajout Olivier Torres
     955    !Config  Key  = choix_bulk
     956    !Config  Desc = choix de la formulation bulk a prendre dans clcdrag au-dessus de l'ocean
     957    !Config  Def  = 0
     958    !Config         0 -> originale (lmdz/Louis 79)
     959    !Config         1 -> COARE
     960    !Config         2 -> CORE-"pure" (cf. Large)
     961    !Config         3 -> CORE-"mixte" (avec z_0 et C_T^N donnees par Smith 88)
     962    choix_bulk_omp = 0
     963    CALL getin('choix_bulk',choix_bulk_omp)
     964
     965    !Config  Key  = nit_bulk
     966    !Config  Desc = choix du nombre d'it de pt fixe dans la bulk
     967    !Config  Def  = 5
     968    nit_bulk_omp = 1
     969    CALL getin('nit_bulk',nit_bulk_omp)
     970
     971    !Config  Key  = kz0
     972    !Config  Desc = choix de la formulation z0 pour la bulk ECUME
     973    !Config  Def  = 1
     974    !Config         0 -> ARPEGE formulation
     975    !Config         1 -> Smith Formulation
     976    !Config         2 -> Direct computation using the stability functions
     977    kz0_omp = 0
     978    CALL getin('kz0',kz0_omp)
     979
    949980
    950981    !Config  Key  = iflag_con
    … …  
    26922723    read_fco2_land_cor = read_fco2_land_cor_omp
    26932724    var_fco2_land_cor = var_fco2_land_cor_omp
     2725    !rajout Olivier Torres
     2726    kz0=kz0_omp
     2727    choix_bulk = choix_bulk_omp
     2728    nit_bulk = nit_bulk_omp
    26942729
    26952730    ! Test of coherence between type_ocean and version_ocean
    … …  
    30563091    WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol
    30573092    WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol
     3093    !rajout Olivier Torres
     3094    write(lunout,*) 'choix_bulk = ', choix_bulk
     3095    write(lunout,*) 'nit_bulk = ', nit_bulk
     3096    write(lunout,*) 'kz0 = ', kz0
    30583097
    30593098    !$OMP END MASTER

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/pbl_surface_mod.F90

    r4478 r4661  
    349349    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
    350350    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
    351     REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
     351    REAL, DIMENSION(klon),        INTENT(INOUT)     :: rain_f  ! rain fall
    352352    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
    353353    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
    … …  
    15701570        ENDDO
    15711571        CALL cdrag(knon, nsrf, &
    1572             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
     1572            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &
    15731573            yts, yqsurf, yz0m, yz0h, yri0, 0, &
    1574             ycdragm, ycdragh, zri1, pref )
     1574            ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1))
    15751575
    15761576! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
    … …  
    16041604
    16051605            CALL cdrag(knon, nsrf, &
    1606             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
     1606            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,&
    16071607            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
    1608             ycdragm_x, ycdragh_x, zri1_x, pref_x )
     1608            ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) )
    16091609
    16101610! --- special Dice. JYG+MPL 25112013
    … …  
    16311631        ENDDO
    16321632        CALL cdrag(knon, nsrf, &
    1633             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
     1633            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&
    16341634            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
    1635             ycdragm_w, ycdragh_w, zri1_w, pref_w )
     1635            ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) )
    16361636!
    16371637        zgeo1(:) = wake_s(:)*zgeo1_w(:) + (1.-wake_s(:))*zgeo1_x(:)
    … …  
    20372037               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
    20382038               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
    2039                yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     2039               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
    20402040          ENDIF
    20412041         
    … …  
    29972997            uzon, vmer, tair1, qair1, zgeo1, &
    29982998            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    2999             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     2999            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
    30003000        ENDIF
    30013001       ELSE  !(iflag_split .eq.0)
    … …  
    30153015            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
    30163016            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    3017             yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
     3017            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
    30183018        CALL stdlevvar(klon, knon, nsrf, zxli, &
    30193019            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
    30203020            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    3021             yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
     3021            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
    30223022        ENDIF
    30233023!!!

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/physiq_mod.F90

    r4489 r4661  
    12121212    REAL qql1(klon),qql2(klon),corrqql
    12131213
     1214
     1215    REAL, dimension(klon) :: pr_et
     1216    REAL, dimension(klon) :: w_et, jlr_g_c, jlr_g_s
     1217
    12141218    REAL pi
    12151219
    … …  
    26672671       ELSE IF (iflag_gusts==2) THEN
    26682672          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
     2673       !!!! modif olivier torres
     2674       ELSE IF (iflag_gusts==3) THEN
     2675          w_et=wstar(1,3)
     2676          jlr_g_s=(0.65*w_et)**2
     2677          pr_et=rain_con*8640
     2678          jlr_g_c = (((19.8*(pr_et(1:klon)**2))/(1.5+pr_et(1:klon)+pr_et(1:klon)**2))**(0.4))**2
     2679          gustiness(1:klon)=jlr_g_c+jlr_g_s
     2680!!       write(*,*) "rain ",pr_et
     2681!!       write(*,*) "jlr_g_c",jlr_g_c
     2682!!       write(*,*) "wstar",wstar(1,3)
     2683!!       write(*,*) "jlr_g_s",jlr_g_s
    26692684          ! ELSE IF (iflag_gusts==2) THEN
    26702685          !    do i = 1, klon

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/screenc_mod.F90

    r4478 r4661  
    1818                         ts, qsurf, z0m, z0h, psol, &
    1919                         ustar, testar, qstar, okri, ri1, &
    20                          pref, delu, delte, delq)
     20                         pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
    2121      IMPLICIT NONE
    2222!-----------------------------------------------------------------------
    … …  
    6060      REAL, dimension(klon), intent(in) :: speed, temp, q_zref
    6161      REAL, intent(in) :: zref
    62       REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol
    63       REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1         
     62      REAL, dimension(klon), intent(IN) :: ts
     63      REAL, dimension(klon), intent(in) :: qsurf, psol
     64      REAL, dimension(klon), intent(inout):: z0m, z0h
     65      REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1   
     66
     67      REAL, dimension(klon), intent(inout) :: s_pblh   
     68      REAL, dimension(klon), intent(in) :: prain   
     69      REAL, dimension(klon), intent(in) :: tsol   
     70      REAL, DIMENSION(klon), INTENT(IN)    :: pat1 !pression premier lev     
    6471!
    6572      REAL, dimension(klon), intent(out) :: pref, delu, delte, delq
    … …  
    8895      CALL cdrag (knon, nsrf, &
    8996                    speed, temp, q_zref, gref, &
    90                     psol, ts, qsurf, z0m, z0h, &
     97                    psol, s_pblh, ts, qsurf, z0m, z0h, &
    9198                    zri_zero,0,                &
    92                     cdram, cdrah, zri1, pref)
     99                    cdram, cdrah, zri1, pref, prain, tsol, pat1)
    93100      DO i = 1, knon
    94101        IF(ok_prescr_ust) THEN
    … …  
    114121                         cdrm, cdrh,  okri, &
    115122                         ri1, iri1, &
    116                          pref, delm, delh, zri1)
     123                         pref, delm, delh, zri1, s_pblh, prain, tsol, pat1)
    117124      IMPLICIT NONE
    118125!-----------------------------------------------------------------------
    … …  
    156163      REAL, dimension(klon), intent(in) :: speed, temp, q_zref
    157164      REAL, intent(in) :: zref
    158       REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol
     165      REAL, dimension(klon), intent(in) :: ts, qsurf, psol
     166      REAL, dimension(klon), intent(inout) :: z0m, z0h
    159167      REAL, dimension(klon), intent(in) :: cdrm, cdrh, ri1         
     168      REAL, dimension(klon), intent(inout) :: s_pblh   
     169      REAL, dimension(klon), intent(in) :: prain   
     170      REAL, dimension(klon), intent(in) :: tsol   
     171      REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev     
    160172      INTEGER, INTENT(IN)  :: iri1 ! Richardson de la 1ere couche
    161173!
    … …  
    180192      CALL cdrag(knon, nsrf, &
    181193                    speed, temp, q_zref, gref, &
    182                     psol, ts, qsurf, z0m, z0h, &
     194                    psol, s_pblh, ts, qsurf, z0m, z0h, &
    183195                    ri1, iri1, &
    184                     cdram, cdrah, zri1, pref)
     196                    cdram, cdrah, zri1, pref, prain, tsol, pat1)
    185197      DO i = 1, knon
    186198        delm(i) = sqrt(cdrm(i))/sqrt(cdram(i))

  • LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/stdlevvar_mod.F90

    r4478 r4661  
    1919                           u1, v1, t1, q1, z1, &
    2020                           ts1, qsurf, z0m, z0h, psol, pat1, &
    21                            t_2m, q_2m, t_10m, q_10m, u_10m, ustar)
     21                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol)
    2222      IMPLICIT NONE
    2323!-------------------------------------------------------------------------
    … …  
    6161      LOGICAL, intent(in) :: zxli
    6262      REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
    63       REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h
     63      REAL, dimension(klon), intent(in) :: qsurf
     64      REAL, dimension(klon), intent(inout) :: z0m, z0h
    6465      REAL, dimension(klon), intent(in) :: psol, pat1
    6566!
    6667      REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar
    6768      REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
     69      REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh
     70      REAL, DIMENSION(klon), INTENT(IN) :: prain
     71      REAL, DIMENSION(klon), INTENT(IN) :: tsol
    6872!-------------------------------------------------------------------------
    6973      include "flux_arp.h"
    … …  
    120124      CALL cdrag(knon, nsrf, &
    121125 &                   speed, t1, q1, z1, &
    122  &                   psol, ts1, qsurf, z0m, z0h, &
     126 &                   psol, s_pblh, ts1, qsurf, z0m, z0h, &
    123127 &                   zri_zero, 0, &
    124  &                   cdram, cdrah, zri1, pref)
     128 &                   cdram, cdrah, zri1, pref, prain, tsol, pat1)
    125129
    126130! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
    … …  
    178182 &                   ts1, qsurf, z0m, z0h, psol, &           
    179183 &                   ustar, testar, qstar, okri, ri1, &
    180  &                   pref, delu, delte, delq)
     184 &                   pref, delu, delte, delq, s_pblh ,prain, tsol, pat1)
    181185!
    182186        DO i = 1, knon
    … …  
    280284 &                   ts1, qsurf, z0m, z0h, psol, &
    281285 &                   ustar, testar, qstar, okri, ri1, &
    282  &                   pref, delu, delte, delq)
     286 &                   pref, delu, delte, delq, s_pblh ,prain, tsol, pat1)
    283287!
    284288        DO i = 1, knon
    … …  
    357361      INTEGER, intent(in) :: klon, knon, nsrf
    358362      LOGICAL, intent(in) :: zxli
    359       REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
    360       REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h
     363      REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, ts1, z1
     364      REAL, dimension(klon), intent(inout) :: z0m, z0h
     365      REAL, dimension(klon), intent(in) :: qsurf
    361366      REAL, dimension(klon), intent(in) :: psol, pat1
    362367!
    … …  
    371376      REAL, dimension(klon) :: cdmn2m, cdhn2m, fm2m, fh2m
    372377      REAL, dimension(klon) :: ri2m_new
     378      REAL, DIMENSION(klon) :: s_pblh
     379      REAL, DIMENSION(klon) :: prain
     380      REAL, DIMENSION(klon) :: tsol
    373381!-------------------------------------------------------------------------
    374382      include "flux_arp.h"
    … …  
    444452      CALL cdrag(knon, nsrf, &
    445453 &                   speed, t1, q1, z1, &
    446  &                   psol, ts1, qsurf, z0m, z0h, &
     454 &                   psol, s_pblh, ts1, qsurf, z0m, z0h, &
    447455 &                   zri_zero, 0, &
    448  &                   cdram, cdrah, zri1, pref)
     456 &                   cdram, cdrah, zri1, pref, prain, tsol, pat1)
    449457
    450458!
    … …  
    468476 &                   cdram, cdrah,  okri, &
    469477 &                   ri1, 1, &
    470  &                   pref_new, delm_new, delh_new, ri2m)
     478 &                   pref_new, delm_new, delh_new, ri2m, &
     479 &                   s_pblh, prain, tsol, pat1      )
    471480!
    472481       DO i = 1, knon
    … …  
    535544 &                   cdram, cdrah,  okri, &
    536545 &                   ri1, 0, &
    537  &                   pref, delm, delh, ri2m)
     546 &                   pref, delm, delh, ri2m, &
     547 &                   s_pblh, prain, tsol, pat1      )
    538548!
    539549        DO i = 1, knon
    … …  
    614624 &                   cdram, cdrah,  okri, &
    615625 &                   ri1, 1, &
    616  &                   pref_new, delm_new, delh_new, ri10m)
     626 &                   pref_new, delm_new, delh_new, ri10m, &
     627 &                   s_pblh, prain, tsol, pat1      )
    617628!
    618629       DO i = 1, knon
    … …  
    671682 &                   cdram, cdrah,  okri, &
    672683 &                   ri1, 0, &
    673  &                   pref, delm, delh, ri10m)
     684 &                   pref, delm, delh, ri10m, &
     685 &                   s_pblh, prain, tsol, pat1      )
    674686!
    675687        DO i = 1, knon
Note: See TracChangeset for help on using the changeset viewer.