Changeset 4722 for LMDZ6


Ignore:
Timestamp:
Oct 9, 2023, 5:33:07 PM (15 months ago)
Author:
Laurent Fairhead
Message:

Modification by O. Torres to the cdrag routines to include different bulk formulae
to calculate cdrag coefficients over ocean as well as an iteration of that
calculation.
The iteration is controlled by flag ok_cdrag_iter which if set to FALSE by default
to converge with previous results.
The choice of bulk formulae is set with the choix_bulk parameter
The number of iterations to run is set with nit_bulk
OT, PB, CD, LF

Location:
LMDZ6/trunk
Files:
8 edited
9 copied

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk

  • LMDZ6/trunk/libf/phylmd/cdrag_mod.F90

    r4687 r4722  
    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
    … …  
    116122  REAL                   MU, CM, CH, B, CMstar, CHstar
    117123  REAL                   PM, PH, BPRIME
    118   REAL                   C
    119124  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
    120125  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
    121   INTEGER                i
     126  INTEGER                i, k
    122127  REAL                   zdu2, ztsolv
    123128  REAL                   ztvd, zscf
    124129  REAL                   zucf, zcr
    125   REAL                   friv, frih
    126130  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
    127131  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
    … …  
    129133  REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
    130134  REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl
     135
     136!------------------ Rajout (OT2018) --------------------
     137!------------------ Rajout pour les appelles BULK (OT) --------------------
     138  REAL, DIMENSION(klon,2) :: rugos_itm
     139  REAL, DIMENSION(klon,2) :: rugos_ith
     140  REAL, PARAMETER         :: tol_it_rugos=1.e-4
     141  REAL, DIMENSION(3)      :: coeffs
     142  LOGICAL                 :: mixte
     143  REAL                    :: z_0m
     144  REAL                    :: z_0h
     145  REAL, DIMENSION(klon)   :: cdmm
     146  REAL, DIMENSION(klon)   :: cdhh
     147
     148!------------------RAJOUT POUR ECUME -------------------
     149
     150  REAL, DIMENSION(klon) :: PQSAT
     151  REAL, DIMENSION(klon) :: PSFTH
     152  REAL, DIMENSION(klon) :: PFSTQ
     153  REAL, DIMENSION(klon) :: PUSTAR
     154  REAL, DIMENSION(klon) :: PCD      ! Drag coefficient for momentum
     155  REAL, DIMENSION(klon) :: PCDN     ! Drag coefficient for momentum
     156  REAL, DIMENSION(klon) :: PCH      ! Drag coefficient for momentum
     157  REAL, DIMENSION(klon) :: PCE      ! Drag coefficient for momentum
     158  REAL, DIMENSION(klon) :: PRI
     159  REAL, DIMENSION(klon) :: PRESA
     160  REAL, DIMENSION(klon) :: PSSS
     161
     162  LOGICAL               :: OPRECIP
     163  LOGICAL               :: OPWEBB
     164  LOGICAL               :: OPERTFLUX
     165  LOGICAL               :: LPRECIP
     166  LOGICAL               :: LPWG
     167
     168
    131169
    132170  LOGICAL, SAVE :: firstcall = .TRUE.
    … …  
    136174  INTEGER, SAVE :: iflag_corr_insta
    137175  !$OMP THREADPRIVATE(iflag_corr_insta)
     176  LOGICAL, SAVE :: ok_cdrag_iter
     177  !$OMP THREADPRIVATE(ok_cdrag_iter)
    138178
    139179!===================================================================c
    … …  
    170210
    171211! Consistent with atke scheme
    172 cn=(1./sqrt(cepsilon))**(2/3)
    173 ri0=2./rpi*(cinf - cn)*ric/cn
    174 ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
     212  cn=(1./sqrt(cepsilon))**(2./3.)
     213  ri0=2./rpi*(cinf - cn)*ric/cn
     214  ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
    175215
    176216
    … …  
    209249! On choisit les fonctions de stabilite utilisees au premier appel
    210250!**************************************************************************
    211   IF (firstcall) THEN
     251 IF (firstcall) THEN
    212252   iflag_corr_sta=2
    213253   iflag_corr_insta=2
     254   ok_cdrag_iter = .FALSE.
    214255 
    215256   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
    216257   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
     258   CALL getin_p('ok_cdrag_iter',ok_cdrag_iter)
    217259
    218260   firstcall = .FALSE.
    219261 ENDIF
     262
     263!------------------ Rajout (OT2018) --------------------
     264!---------    Rajout pour itération sur rugosité     ----------------
     265  rugos_itm(:,2) = z0m
     266  rugos_itm(:,1) = 3.*tol_it_rugos*z0m
     267
     268  rugos_ith(:,2) = z0h !cp nouveau rugos_it
     269  rugos_ith(:,1) = 3.*tol_it_rugos*z0h
     270!--------------------------------------------------------------------
    220271
    221272!xxxxxxxxxxxxxxxxxxxxxxx
    … …  
    227278!***********************
    228279     
    229 
    230280     zdu2 = MAX(CEPDU2, speed(i)**2)
    231281     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
    232                  (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
     282          (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
    233283     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))) &
     284     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
    235285          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
    236      zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     286     
     287!------------------ Rajout (OT2018) --------------------
     288!--------------   ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN    ------------------------     
    237289     IF (iri_in.EQ.1) THEN
    238290       zri(i) = ri_in(i)
    239291     ENDIF
    240292
     293     IF (nsrf == is_oce) THEN
     294       
     295!------------------  Pour Core 2 choix Core Pur ou Core Mixte  --------------------
     296       IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN
     297         IF(choix_bulk .eq. 2) THEN
     298           mixte = .false.
     299         ELSE
     300            mixte = .true.
     301         ENDIF
     302         call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),&
     303             zgeop1(i)/RG, zgeop1(i)/RG,  zgeop1(i)/RG,&
     304             psol(i),nit_bulk,mixte,&
     305             coeffs,z_0m,z_0h)
     306         cdmm(i) = coeffs(1)
     307         cdhh(i) = coeffs(2)
     308         cdm(i)=cdmm(i)
     309         cdh(i)=cdhh(i)
     310         write(*,*) "clc_core cd ch",cdmm(i),cdhh(i)
     311
     312!------------------                 Pour ECUME                 --------------------
     313       ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN
     314         OPRECIP = .false.
     315         OPWEBB  = .false.
     316         OPERTFLUX = .false.
     317         IF (nsrf .eq. is_oce) THEN
     318           PSSS = 0.0
     319         ENDIF
     320         call ini_csts
     321         call ecumev6_flux( z_0m,t1(i),tsurf(i),&
     322             q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,&
     323             psol(i),pat1(i), OPRECIP, OPWEBB,&
     324             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,&
     325             PCH,PCE,PRI,PRESA,prain,&
     326             z_0h,OPERTFLUX,coeffs)
     327
     328         cdmm(i) = coeffs(1)
     329         cdhh(i) = coeffs(2)
     330         cdm(i)=cdmm(i)
     331         cdh(i)=cdhh(i)
     332   
     333         write(*,*) "ecume cd ch",cdmm(i),cdhh(i)
     334
     335!------------------                 Pour COARE CNRM                 --------------------
     336       ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN
     337         LPRECIP = .false.
     338         LPWG    = .false.
     339         call ini_csts
     340         call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i),  &
     341             sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, &
     342             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
     343             PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs)
     344         cdmm(i) = coeffs(1)
     345         cdhh(i) = coeffs(2)
     346         cdm(i)=cdmm(i)
     347         cdh(i)=cdhh(i)
     348         write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i)
     349
     350!------------------                 Pour COARE Maison                 --------------------
     351       ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN
     352         IF ( pblh(i) .eq. 0. ) THEN
     353           pblh(i) = 1500.
     354         ENDIF
     355         write(*,*) "debug size",size(coeffs)
     356         call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),&
     357               t1(i),q1(i),&
     358               zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,&
     359               psol(i), pblh(i),&
     360               nit_bulk,&
     361               coeffs,z_0m,z_0h)
     362         cdmm(i) = coeffs(1)
     363         cdhh(i) = coeffs(3)
     364         cdm(i)=cdmm(i)
     365         cdh(i)=cdhh(i)
     366         write(*,*) "coare cd ch",cdmm(i),cdhh(i)
     367       ELSE
     368!------------------                 Pour La param LMDZ (ocean)              --------------------
     369         zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     370         IF (iri_in.EQ.1) THEN
     371           zri(i) = ri_in(i)
     372         ENDIF
     373       ENDIF
     374     
     375
     376!----------------------- Rajout des itérations --------------
     377       DO  k=1,nit_bulk
    241378
    242379! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
    243380!********************************************************************
    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))))
     381         zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2)))
     382         cdmn(i) = zzzcd*zzzcd
     383         cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2))))
     384
     385! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
     386!*******************************************************
     387         IF (zri(i) .LT. 0.) THEN   
     388           SELECT CASE (iflag_corr_insta)
     389             CASE (1) ! Louis 1979 + Mascart 1995
     390                  MU=LOG(MAX(z0m(i)/z0h(i),0.01))
     391                  CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
     392                  PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
     393                  CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
     394                  PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
     395                  CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     396                     & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
     397                     & * ((zgeop1(i)/(RG*z0h(i)))**PH)
     398                  CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     399                     & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
     400                     & * ((zgeop1(i)/(RG*z0m(i)))**PM)
     401                  FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
     402                  FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
     403             CASE (2) ! Louis 1982
     404                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     405                        *(1.0+zgeop1(i)/(RG*z0m(i)))))
     406                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     407                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     408             CASE (3) ! Laurent Li
     409                  FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     410                  FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     411             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     412                  sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
     413                  prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
     414                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
     415                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     416             CASE default ! Louis 1982
     417                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     418                         *(1.0+zgeop1(i)/(RG*z0m(i)))))
     419                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     420                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     421           END SELECT
     422! Calcul des drags
     423           cdmm(i)=cdmn(i)*FM(i)
     424           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     425! Traitement particulier des cas oceaniques
     426! on applique Miller et al 1992 en l'absence de gustiness
     427           IF (nsrf == is_oce) THEN
     428!            cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     429             IF (iflag_gusts==0) THEN
     430               zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
     431               cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
     432             ENDIF
     433             cdmm(i)=MIN(cdmm(i),cdmmax)
     434             cdhh(i)=MIN(cdhh(i),cdhmax)
     435!             write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i)
     436           END IF
     437         ELSE                         
     438
     439!'''''''''''''''
     440! Cas stables :
     441!'''''''''''''''
     442           zri(i) = MIN(20.,zri(i))
     443           SELECT CASE (iflag_corr_sta)
     444             CASE (1) ! Louis 1979 + Mascart 1995
     445                  FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
     446                  FH(i)=FM(i)
     447             CASE (2) ! Louis 1982
     448                  zscf = SQRT(1.+CD*ABS(zri(i)))
     449                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     450                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     451             CASE (3) ! Laurent Li
     452                  FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
     453                  FH(i)=FM(i)
     454             CASE (4)  ! King 2001
     455                  IF (zri(i) .LT. C2/2.) THEN
     456                    FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
     457                    FH(i)=  FM(i)
     458                  ELSE
     459                    FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
     460                    FH(i)= FM(i)
     461                  ENDIF
     462             CASE (5) ! MO
     463                  if (zri(i) .LT. 1./alpha) then
     464                    FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
     465                    FH(i)=FM(i)
     466                  else
     467                    FM(i)=MAX(1E-7,f_ri_cd_min)
     468                    FH(i)=FM(i)
     469                  endif
     470             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     471                  sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
     472                  prandtl(i) = pr_neut + zri(i) * pr_slope
     473                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
     474                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     475             CASE default ! Louis 1982
     476                  zscf = SQRT(1.+CD*ABS(zri(i)))
     477                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     478                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     479           END SELECT
     480
     481! Calcul des drags
     482           
     483           cdmm(i)=cdmn(i)*FM(i)
     484           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     485           IF (choix_bulk == 0) THEN
     486             cdm(i)=cdmn(i)*FM(i)
     487             cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     488           ENDIF
     489
     490           IF (nsrf.EQ.is_oce) THEN
     491             cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i)
     492             cdmm(i)=MIN(cdm(i),cdmmax)
     493             cdhh(i)=MIN(cdh(i),cdhmax)
     494           ENDIF
     495           IF (ok_cdrag_iter) THEN
     496             rugos_itm(i,1) = rugos_itm(i,2)
     497             rugos_ith(i,1) = rugos_ith(i,2)
     498             rugos_itm(i,2) =  0.018*cdmm(i) * (speed(i))/RG  &
     499                              +  0.11*14e-6 / SQRT(cdmm(i) * zdu2)
     500
     501!---------- Version SEPARATION DES Z0 ----------------------
     502             IF (iflag_z0_oce==0) THEN
     503               rugos_ith(i,2) = rugos_itm(i,2)
     504             ELSE IF (iflag_z0_oce==1) THEN
     505               rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2)
     506             ENDIF
     507           ENDIF
     508         ENDIF
     509         IF (ok_cdrag_iter) THEN
     510           rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2))
     511           rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2))
     512         ENDIF
     513       ENDDO
     514       IF (nsrf.EQ.is_oce) THEN
     515         cdm(i)=MIN(cdmm(i),cdmmax)
     516         cdh(i)=MIN(cdhh(i),cdhmax)
     517       ENDIF
     518       z0m = rugos_itm(:,2)
     519       z0h = rugos_ith(:,2)
     520     ELSE ! (nsrf == is_oce)
     521       zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
     522       IF (iri_in.EQ.1) THEN
     523         zri(i) = ri_in(i)
     524       ENDIF
     525
     526! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
     527!********************************************************************
     528       zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
     529       cdmn(i) = zzzcd*zzzcd
     530       cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
    248531
    249532
    250533! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
    251534!*******************************************************
    252 
    253 
    254 
    255535!''''''''''''''
    256536! Cas instables
    257537!''''''''''''''
    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 
     538       IF (zri(i) .LT. 0.) THEN   
     539         SELECT CASE (iflag_corr_insta)
     540           CASE (1) ! Louis 1979 + Mascart 1995
     541                MU=LOG(MAX(z0m(i)/z0h(i),0.01))
     542                CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
     543                PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
     544                CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
     545                PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
     546                CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     547                   & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
     548                   & * ((zgeop1(i)/(RG*z0h(i)))**PH)
     549                CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
     550                   & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
     551                   & * ((zgeop1(i)/(RG*z0m(i)))**PM)
     552                FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
     553                FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
     554           CASE (2) ! Louis 1982
     555                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     556                       *(1.0+zgeop1(i)/(RG*z0m(i)))))
     557                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     558                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
     559           CASE (3) ! Laurent Li
     560                FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     561                FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
     562           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     563                 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
     564                 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
     565                 FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
     566                 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     567           CASE default ! Louis 1982
     568                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
     569                      *(1.0+zgeop1(i)/(RG*z0m(i)))))
     570                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
     571                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
    313572         END SELECT
    314 
    315 
    316 
    317573! 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 
     574         cdm(i)=cdmn(i)*FM(i)
     575         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     576       ELSE                         
    345577!'''''''''''''''
    346578! Cas stables :
    347579!'''''''''''''''
    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 
     580         zri(i) = MIN(20.,zri(i))
     581         SELECT CASE (iflag_corr_sta)
     582           CASE (1) ! Louis 1979 + Mascart 1995
     583                FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
     584                FH(i)=FM(i)
     585           CASE (2) ! Louis 1982
     586                zscf = SQRT(1.+CD*ABS(zri(i)))
     587                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     588                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     589           CASE (3) ! Laurent Li
     590                FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
     591                FH(i)=FM(i)
     592           CASE (4)  ! King 2001
     593                if (zri(i) .LT. C2/2.) then
     594                  FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
     595                  FH(i)=  FM(i)
     596                else
     597                  FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
     598                  FH(i)= FM(i)
     599                endif
     600           CASE (5) ! MO
     601                if (zri(i) .LT. 1./alpha) then
     602                  FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
     603                  FH(i)=FM(i)
     604                else
     605                  FM(i)=MAX(1E-7,f_ri_cd_min)
     606                  FH(i)=FM(i)
     607                endif
     608           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
     609                sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
     610                prandtl(i) = pr_neut + zri(i) * pr_slope
     611                FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
     612                FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
     613           CASE default ! Louis 1982
     614                zscf = SQRT(1.+CD*ABS(zri(i)))
     615                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
     616                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
     617         END SELECT
    415618! 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
     619         cdm(i)=cdmn(i)*FM(i)
     620         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
     621       ENDIF
     622     ENDIF ! fin du if (nsrf == is_oce)
    433623  END DO   !  Fin de la boucle sur l'horizontal
    434 !xxxxxxxxxxx
    435 ! ================================================================= c
    436624     
    437625END SUBROUTINE cdrag

  • LMDZ6/trunk/libf/phylmd/clesphys.h

    r4677 r4722  
    4343!IM seuils cdrm, cdrh
    4444       REAL cdmmax, cdhmax
     45!IM pour les params différentes Olivier Torres
     46       INTEGER choix_bulk, nit_bulk, kz0
    4547!IM param. stabilite s/ terres et en dehors
    4648       REAL ksta, ksta_ter, f_ri_cd_min
    … …  
    111113       COMMON/clesphys/                                                 &
    112114! REAL FIRST
     115! rajout choix_bulk et nit_bulk kz0 par Olivier Torres
    113116     &       co2_ppm, solaire                                           &
    114117     &     , RCO2, RCH4, RN2O, RCFC11, RCFC12                           &
    … …  
    138141     &     , ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad                &
    139142     &     , iflag_con, nbapp_cv, nbapp_wk                              &
     143     &     , choix_bulk, nit_bulk, kz0                                  &
    140144     &     , iflag_ener_conserv                                         &
    141145     &     , ok_suntime_rrtm                                            &

  • LMDZ6/trunk/libf/phylmd/conf_phys_m.F90

    r4707 r4722  
    243243    LOGICAL, SAVE :: ok_new_lscp_omp
    244244    LOGICAL, SAVE :: ok_icefra_lscp_omp
     245    !rajout de choix_bulk et nit_bulk par Olivier Torres
     246    INTEGER,SAVE  :: choix_bulk_omp
     247    INTEGER,SAVE  :: nit_bulk_omp
     248    INTEGER,SAVE  :: kz0_omp
    245249    LOGICAL, SAVE :: ok_bs_omp, ok_rad_bs_omp
    246250
    … …  
    936940    nbapp_rad_omp = 12
    937941    CALL getin('nbapp_rad',nbapp_rad_omp)
     942
     943    !rajout Olivier Torres
     944    !Config  Key  = choix_bulk
     945    !Config  Desc = choix de la formulation bulk a prendre dans clcdrag au-dessus de l'ocean
     946    !Config  Def  = 0
     947    !Config         0 -> originale (lmdz/Louis 79)
     948    !Config         1 -> COARE
     949    !Config         2 -> CORE-"pure" (cf. Large)
     950    !Config         3 -> CORE-"mixte" (avec z_0 et C_T^N donnees par Smith 88)
     951    choix_bulk_omp = 0
     952    CALL getin('choix_bulk',choix_bulk_omp)
     953
     954    !Config  Key  = nit_bulk
     955    !Config  Desc = choix du nombre d'it de pt fixe dans la bulk
     956    !Config  Def  = 5
     957    nit_bulk_omp = 1
     958    CALL getin('nit_bulk',nit_bulk_omp)
     959
     960    !Config  Key  = kz0
     961    !Config  Desc = choix de la formulation z0 pour la bulk ECUME
     962    !Config  Def  = 1
     963    !Config         0 -> ARPEGE formulation
     964    !Config         1 -> Smith Formulation
     965    !Config         2 -> Direct computation using the stability functions
     966    kz0_omp = 0
     967    CALL getin('kz0',kz0_omp)
     968
    938969
    939970    !Config  Key  = iflag_con
    … …  
    24962527    var_fco2_land_cor = var_fco2_land_cor_omp
    24972528    dms_cycle_cpl = dms_cycle_cpl_omp
     2529    !rajout Olivier Torres
     2530    kz0=kz0_omp
     2531    choix_bulk = choix_bulk_omp
     2532    nit_bulk = nit_bulk_omp
    24982533
    24992534    ! Test of coherence between type_ocean and version_ocean
    … …  
    28472882    WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol
    28482883    WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol
     2884    !rajout Olivier Torres
     2885    write(lunout,*) 'choix_bulk = ', choix_bulk
     2886    write(lunout,*) 'nit_bulk = ', nit_bulk
     2887    write(lunout,*) 'kz0 = ', kz0
    28492888
    28502889    !$OMP END MASTER

  • LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

    r4687 r4722  
    348348    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
    349349    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
    350     REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
     350    REAL, DIMENSION(klon),        INTENT(INOUT)     :: rain_f  ! rain fall
    351351    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
    352352    REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
    … …  
    10711071!albedo SB <<<
    10721072    yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
    1073     ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0   
     1073    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yz0h_oupas = 0.0 ; yu1 = 0.0   
    10741074    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
    10751075    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
    … …  
    15981598        ENDDO
    15991599        CALL cdrag(knon, nsrf, &
    1600             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
     1600            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &
    16011601            yts, yqsurf, yz0m, yz0h, yri0, 0, &
    1602             ycdragm, ycdragh, zri1, pref )
     1602            ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1))
    16031603
    16041604! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
    … …  
    16321632
    16331633            CALL cdrag(knon, nsrf, &
    1634             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
     1634            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,&
    16351635            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
    1636             ycdragm_x, ycdragh_x, zri1_x, pref_x )
     1636            ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) )
    16371637
    16381638! --- special Dice. JYG+MPL 25112013
    … …  
    16591659        ENDDO
    16601660        CALL cdrag(knon, nsrf, &
    1661             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
     1661            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&
    16621662            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
    1663             ycdragm_w, ycdragh_w, zri1_w, pref_w )
     1663            ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) )
    16641664!
    16651665!!!bug !!        zgeo1(:) = wake_s(:)*zgeo1_w(:) + (1.-wake_s(:))*zgeo1_x(:)
    … …  
    20862086               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
    20872087               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
    2088                yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     2088               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
    20892089          ENDIF
    20902090         
    … …  
    30723072       IF (iflag_split .eq.0) THEN
    30733073        IF (iflag_new_t2mq2m==1) THEN
    3074          CALL stdlevvarn(klon, knon, nsrf, zxli, &
     3074           CALL stdlevvarn(klon, knon, nsrf, zxli, &
    30753075            uzon, vmer, tair1, qair1, zgeo1, &
    30763076            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    … …  
    30813081            uzon, vmer, tair1, qair1, zgeo1, &
    30823082            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    3083             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
     3083            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
    30843084        ENDIF
    30853085       ELSE  !(iflag_split .eq.0)
    … …  
    30993099            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
    31003100            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    3101             yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
     3101            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
    31023102        CALL stdlevvar(klon, knon, nsrf, zxli, &
    31033103            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
    31043104            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
    3105             yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
     3105            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
    31063106        ENDIF
    31073107!!!

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4715 r4722  
    12591259    REAL, dimension(klon,klev) :: t_env,q_env
    12601260
     1261    REAL, dimension(klon) :: pr_et
     1262    REAL, dimension(klon) :: w_et, jlr_g_c, jlr_g_s
     1263
    12611264    REAL pi
    12621265    INTEGER ieru
    … …  
    28042807       ELSE IF (iflag_gusts==2) THEN
    28052808          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
     2809       !!!! modif olivier torres
     2810       ELSE IF (iflag_gusts==3) THEN
     2811          w_et=wstar(1,3)
     2812          jlr_g_s=(0.65*w_et)**2
     2813          pr_et=rain_con*8640
     2814          jlr_g_c = (((19.8*(pr_et(1:klon)**2))/(1.5+pr_et(1:klon)+pr_et(1:klon)**2))**(0.4))**2
     2815          gustiness(1:klon)=jlr_g_c+jlr_g_s
     2816!!       write(*,*) "rain ",pr_et
     2817!!       write(*,*) "jlr_g_c",jlr_g_c
     2818!!       write(*,*) "wstar",wstar(1,3)
     2819!!       write(*,*) "jlr_g_s",jlr_g_s
    28062820          ! ELSE IF (iflag_gusts==2) THEN
    28072821          !    do i = 1, klon

  • LMDZ6/trunk/libf/phylmd/screenc_mod.F90

    r4478 r4722  
    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/trunk/libf/phylmd/stdlevvar_mod.F90

    r4478 r4722  
    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.