Index: LMDZ6/trunk/libf/phylmd/cdrag_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/cdrag_mod.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/cdrag_mod.F90 (revision 4722) @@ -15,7 +15,7 @@ SUBROUTINE cdrag(knon, nsrf, & speed, t1, q1, zgeop1, & - psol, tsurf, qsurf, z0m, z0h, & + psol, pblh, tsurf, qsurf, z0m, z0h, & ri_in, iri_in, & - cdm, cdh, zri, pref) + cdm, cdh, zri, pref, prain, tsol , pat1) USE dimphy @@ -81,14 +81,20 @@ INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface - REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele - REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele - REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) - REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) - REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) + REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele + REAL, DIMENSION(klon), INTENT(IN) :: zgeop1 ! geopotentiel au 1er niveau du modele + REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) + REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) + REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h ! Rugosity at surface (m) REAL, DIMENSION(klon), INTENT(IN) :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var. INTEGER, INTENT(IN) :: iri_in! iflag to activate cdrag calculation using ri1 - REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) - REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) - REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol + REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) + REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) + REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol + +!------------------ Rajout pour COARE (OT2018) -------------------- + REAL, DIMENSION(klon), INTENT(INOUT) :: pblh !hauteur de CL + REAL, DIMENSION(klon), INTENT(IN) :: prain !rapport au precipitation + REAL, DIMENSION(klon), INTENT(IN) :: tsol !SST imposé sur la surface oceanique + REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev @@ -99,5 +105,5 @@ REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for heat flux REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number - REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG + REAL, DIMENSION(klon), INTENT(INOUT) :: pref ! Pression au niveau zgeop/RG ! Variables Locales @@ -116,12 +122,10 @@ REAL MU, CM, CH, B, CMstar, CHstar REAL PM, PH, BPRIME - REAL C INTEGER ng_q1 ! Number of grids that q1 < 0.0 INTEGER ng_qsurf ! Number of grids that qsurf < 0.0 - INTEGER i + INTEGER i, k REAL zdu2, ztsolv REAL ztvd, zscf REAL zucf, zcr - REAL friv, frih REAL, DIMENSION(klon) :: FM, FH ! stability functions REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions @@ -129,4 +133,38 @@ REAL, DIMENSION(klon) :: sm, prandtl ! Stability function from atke turbulence scheme REAL ri0, ri1, cn ! to have dimensionless term in sm and prandtl + +!------------------ Rajout (OT2018) -------------------- +!------------------ Rajout pour les appelles BULK (OT) -------------------- + REAL, DIMENSION(klon,2) :: rugos_itm + REAL, DIMENSION(klon,2) :: rugos_ith + REAL, PARAMETER :: tol_it_rugos=1.e-4 + REAL, DIMENSION(3) :: coeffs + LOGICAL :: mixte + REAL :: z_0m + REAL :: z_0h + REAL, DIMENSION(klon) :: cdmm + REAL, DIMENSION(klon) :: cdhh + +!------------------RAJOUT POUR ECUME ------------------- + + REAL, DIMENSION(klon) :: PQSAT + REAL, DIMENSION(klon) :: PSFTH + REAL, DIMENSION(klon) :: PFSTQ + REAL, DIMENSION(klon) :: PUSTAR + REAL, DIMENSION(klon) :: PCD ! Drag coefficient for momentum + REAL, DIMENSION(klon) :: PCDN ! Drag coefficient for momentum + REAL, DIMENSION(klon) :: PCH ! Drag coefficient for momentum + REAL, DIMENSION(klon) :: PCE ! Drag coefficient for momentum + REAL, DIMENSION(klon) :: PRI + REAL, DIMENSION(klon) :: PRESA + REAL, DIMENSION(klon) :: PSSS + + LOGICAL :: OPRECIP + LOGICAL :: OPWEBB + LOGICAL :: OPERTFLUX + LOGICAL :: LPRECIP + LOGICAL :: LPWG + + LOGICAL, SAVE :: firstcall = .TRUE. @@ -136,4 +174,6 @@ INTEGER, SAVE :: iflag_corr_insta !$OMP THREADPRIVATE(iflag_corr_insta) + LOGICAL, SAVE :: ok_cdrag_iter + !$OMP THREADPRIVATE(ok_cdrag_iter) !===================================================================c @@ -170,7 +210,7 @@ ! Consistent with atke scheme -cn=(1./sqrt(cepsilon))**(2/3) -ri0=2./rpi*(cinf - cn)*ric/cn -ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope + cn=(1./sqrt(cepsilon))**(2./3.) + ri0=2./rpi*(cinf - cn)*ric/cn + ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope @@ -209,13 +249,24 @@ ! On choisit les fonctions de stabilite utilisees au premier appel !************************************************************************** - IF (firstcall) THEN + IF (firstcall) THEN iflag_corr_sta=2 iflag_corr_insta=2 + ok_cdrag_iter = .FALSE. CALL getin_p('iflag_corr_sta',iflag_corr_sta) CALL getin_p('iflag_corr_insta',iflag_corr_insta) + CALL getin_p('ok_cdrag_iter',ok_cdrag_iter) firstcall = .FALSE. ENDIF + +!------------------ Rajout (OT2018) -------------------- +!--------- Rajout pour itération sur rugosité ---------------- + rugos_itm(:,2) = z0m + rugos_itm(:,1) = 3.*tol_it_rugos*z0m + + rugos_ith(:,2) = z0h !cp nouveau rugos_it + rugos_ith(:,1) = 3.*tol_it_rugos*z0h +!-------------------------------------------------------------------- !xxxxxxxxxxxxxxxxxxxxxxx @@ -227,211 +278,348 @@ !*********************** - zdu2 = MAX(CEPDU2, speed(i)**2) pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & - (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero + (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero - ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) & + ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero - zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) + +!------------------ Rajout (OT2018) -------------------- +!-------------- ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN ------------------------ IF (iri_in.EQ.1) THEN zri(i) = ri_in(i) ENDIF + IF (nsrf == is_oce) THEN + +!------------------ Pour Core 2 choix Core Pur ou Core Mixte -------------------- + IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN + IF(choix_bulk .eq. 2) THEN + mixte = .false. + ELSE + mixte = .true. + ENDIF + call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),& + zgeop1(i)/RG, zgeop1(i)/RG, zgeop1(i)/RG,& + psol(i),nit_bulk,mixte,& + coeffs,z_0m,z_0h) + cdmm(i) = coeffs(1) + cdhh(i) = coeffs(2) + cdm(i)=cdmm(i) + cdh(i)=cdhh(i) + write(*,*) "clc_core cd ch",cdmm(i),cdhh(i) + +!------------------ Pour ECUME -------------------- + ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN + OPRECIP = .false. + OPWEBB = .false. + OPERTFLUX = .false. + IF (nsrf .eq. is_oce) THEN + PSSS = 0.0 + ENDIF + call ini_csts + call ecumev6_flux( z_0m,t1(i),tsurf(i),& + q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,& + psol(i),pat1(i), OPRECIP, OPWEBB,& + PSFTH,PFSTQ,PUSTAR,PCD,PCDN,& + PCH,PCE,PRI,PRESA,prain,& + z_0h,OPERTFLUX,coeffs) + + cdmm(i) = coeffs(1) + cdhh(i) = coeffs(2) + cdm(i)=cdmm(i) + cdh(i)=cdhh(i) + + write(*,*) "ecume cd ch",cdmm(i),cdhh(i) + +!------------------ Pour COARE CNRM -------------------- + ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN + LPRECIP = .false. + LPWG = .false. + call ini_csts + call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i), & + sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, & + PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, & + PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs) + cdmm(i) = coeffs(1) + cdhh(i) = coeffs(2) + cdm(i)=cdmm(i) + cdh(i)=cdhh(i) + write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i) + +!------------------ Pour COARE Maison -------------------- + ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN + IF ( pblh(i) .eq. 0. ) THEN + pblh(i) = 1500. + ENDIF + write(*,*) "debug size",size(coeffs) + call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),& + t1(i),q1(i),& + zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,& + psol(i), pblh(i),& + nit_bulk,& + coeffs,z_0m,z_0h) + cdmm(i) = coeffs(1) + cdhh(i) = coeffs(3) + cdm(i)=cdmm(i) + cdh(i)=cdhh(i) + write(*,*) "coare cd ch",cdmm(i),cdhh(i) + ELSE +!------------------ Pour La param LMDZ (ocean) -------------------- + zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) + IF (iri_in.EQ.1) THEN + zri(i) = ri_in(i) + ENDIF + ENDIF + + +!----------------------- Rajout des itérations -------------- + DO k=1,nit_bulk ! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): !******************************************************************** - - zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) - cdmn(i) = zzzcd*zzzcd - cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) + zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2))) + cdmn(i) = zzzcd*zzzcd + cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2)))) + +! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi : +!******************************************************* + IF (zri(i) .LT. 0.) THEN + SELECT CASE (iflag_corr_insta) + CASE (1) ! Louis 1979 + Mascart 1995 + MU=LOG(MAX(z0m(i)/z0h(i),0.01)) + CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) + PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) + CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) + PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) + CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & + & * ((zgeop1(i)/(RG*z0h(i)))**PH) + CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * ((zgeop1(i)/(RG*z0m(i)))**PM) + FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) + FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) + CASE (2) ! Louis 1982 + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) + CASE (3) ! Laurent Li + FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) + sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn + prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut + FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min) + FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) + CASE default ! Louis 1982 + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) + END SELECT +! Calcul des drags + cdmm(i)=cdmn(i)*FM(i) + cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i) +! Traitement particulier des cas oceaniques +! on applique Miller et al 1992 en l'absence de gustiness + IF (nsrf == is_oce) THEN +! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) + IF (iflag_gusts==0) THEN + zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) + cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) + ENDIF + cdmm(i)=MIN(cdmm(i),cdmmax) + cdhh(i)=MIN(cdhh(i),cdhmax) +! write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i) + END IF + ELSE + +!''''''''''''''' +! Cas stables : +!''''''''''''''' + zri(i) = MIN(20.,zri(i)) + SELECT CASE (iflag_corr_sta) + CASE (1) ! Louis 1979 + Mascart 1995 + FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) + FH(i)=FM(i) + CASE (2) ! Louis 1982 + zscf = SQRT(1.+CD*ABS(zri(i))) + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + CASE (3) ! Laurent Li + FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) + FH(i)=FM(i) + CASE (4) ! King 2001 + IF (zri(i) .LT. C2/2.) THEN + FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) + FH(i)= FM(i) + ELSE + FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) + FH(i)= FM(i) + ENDIF + CASE (5) ! MO + if (zri(i) .LT. 1./alpha) then + FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) + FH(i)=FM(i) + else + FM(i)=MAX(1E-7,f_ri_cd_min) + FH(i)=FM(i) + endif + CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) + sm(i) = MAX(0.,cn*(1.-zri(i)/ric)) + prandtl(i) = pr_neut + zri(i) * pr_slope + FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min) + FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) + CASE default ! Louis 1982 + zscf = SQRT(1.+CD*ABS(zri(i))) + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + END SELECT + +! Calcul des drags + + cdmm(i)=cdmn(i)*FM(i) + cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i) + IF (choix_bulk == 0) THEN + cdm(i)=cdmn(i)*FM(i) + cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) + ENDIF + + IF (nsrf.EQ.is_oce) THEN + cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i) + cdmm(i)=MIN(cdm(i),cdmmax) + cdhh(i)=MIN(cdh(i),cdhmax) + ENDIF + IF (ok_cdrag_iter) THEN + rugos_itm(i,1) = rugos_itm(i,2) + rugos_ith(i,1) = rugos_ith(i,2) + rugos_itm(i,2) = 0.018*cdmm(i) * (speed(i))/RG & + + 0.11*14e-6 / SQRT(cdmm(i) * zdu2) + +!---------- Version SEPARATION DES Z0 ---------------------- + IF (iflag_z0_oce==0) THEN + rugos_ith(i,2) = rugos_itm(i,2) + ELSE IF (iflag_z0_oce==1) THEN + rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2) + ENDIF + ENDIF + ENDIF + IF (ok_cdrag_iter) THEN + rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2)) + rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2)) + ENDIF + ENDDO + IF (nsrf.EQ.is_oce) THEN + cdm(i)=MIN(cdmm(i),cdmmax) + cdh(i)=MIN(cdhh(i),cdhmax) + ENDIF + z0m = rugos_itm(:,2) + z0h = rugos_ith(:,2) + ELSE ! (nsrf == is_oce) + zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) + IF (iri_in.EQ.1) THEN + zri(i) = ri_in(i) + ENDIF + +! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): +!******************************************************************** + zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) + cdmn(i) = zzzcd*zzzcd + cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) ! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi : !******************************************************* - - - !'''''''''''''' ! Cas instables !'''''''''''''' - - IF (zri(i) .LT. 0.) THEN - - - SELECT CASE (iflag_corr_insta) - - CASE (1) ! Louis 1979 + Mascart 1995 - - MU=LOG(MAX(z0m(i)/z0h(i),0.01)) - CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) - PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) - CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) - PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) - CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & - & * ((zgeop1(i)/(RG*z0h(i)))**PH) - CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & * ((zgeop1(i)/(RG*z0m(i)))**PM) - - - - - FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) - FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) - - CASE (2) ! Louis 1982 - - zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & - *(1.0+zgeop1(i)/(RG*z0m(i))))) - FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) - FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) - - - CASE (3) ! Laurent Li - - - FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) - FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) - - CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) - - sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn - prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut - FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min) - FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) - - CASE default ! Louis 1982 - - zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & - *(1.0+zgeop1(i)/(RG*z0m(i))))) - FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) - FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) - - + IF (zri(i) .LT. 0.) THEN + SELECT CASE (iflag_corr_insta) + CASE (1) ! Louis 1979 + Mascart 1995 + MU=LOG(MAX(z0m(i)/z0h(i),0.01)) + CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) + PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) + CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) + PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) + CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & + & * ((zgeop1(i)/(RG*z0h(i)))**PH) + CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * ((zgeop1(i)/(RG*z0m(i)))**PM) + FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) + FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) + CASE (2) ! Louis 1982 + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) + CASE (3) ! Laurent Li + FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) + sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn + prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut + FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min) + FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) + CASE default ! Louis 1982 + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) END SELECT - - - ! Calcul des drags - - - cdm(i)=cdmn(i)*FM(i) - cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) - - -! Traitement particulier des cas oceaniques -! on applique Miller et al 1992 en l'absence de gustiness - - IF (nsrf == is_oce) THEN -! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) - - IF(iflag_gusts==0) THEN - zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) - cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) - ENDIF - - - cdm(i)=MIN(cdm(i),cdmmax) - cdh(i)=MIN(cdh(i),cdhmax) - - END IF - - - - ELSE - + cdm(i)=cdmn(i)*FM(i) + cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) + ELSE !''''''''''''''' ! Cas stables : !''''''''''''''' - zri(i) = MIN(20.,zri(i)) - - SELECT CASE (iflag_corr_sta) - - CASE (1) ! Louis 1979 + Mascart 1995 - - FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) - FH(i)=FM(i) - - - CASE (2) ! Louis 1982 - - zscf = SQRT(1.+CD*ABS(zri(i))) - FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) - FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) - - - CASE (3) ! Laurent Li - - FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) - FH(i)=FM(i) - - - CASE (4) ! King 2001 - - if (zri(i) .LT. C2/2.) then - FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) - FH(i)= FM(i) - - - else - FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) - FH(i)= FM(i) - endif - - - CASE (5) ! MO - - if (zri(i) .LT. 1./alpha) then - - FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) - FH(i)=FM(i) - - else - - - FM(i)=MAX(1E-7,f_ri_cd_min) - FH(i)=FM(i) - - endif - - CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) - sm(i) = MAX(0.,cn*(1.-zri(i)/ric)) - prandtl(i) = pr_neut + zri(i) * pr_slope - FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min) - FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) - - CASE default ! Louis 1982 - - zscf = SQRT(1.+CD*ABS(zri(i))) - FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) - FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) - - - - END SELECT - + zri(i) = MIN(20.,zri(i)) + SELECT CASE (iflag_corr_sta) + CASE (1) ! Louis 1979 + Mascart 1995 + FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) + FH(i)=FM(i) + CASE (2) ! Louis 1982 + zscf = SQRT(1.+CD*ABS(zri(i))) + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + CASE (3) ! Laurent Li + FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) + FH(i)=FM(i) + CASE (4) ! King 2001 + if (zri(i) .LT. C2/2.) then + FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) + FH(i)= FM(i) + else + FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) + FH(i)= FM(i) + endif + CASE (5) ! MO + if (zri(i) .LT. 1./alpha) then + FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) + FH(i)=FM(i) + else + FM(i)=MAX(1E-7,f_ri_cd_min) + FH(i)=FM(i) + endif + CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) + sm(i) = MAX(0.,cn*(1.-zri(i)/ric)) + prandtl(i) = pr_neut + zri(i) * pr_slope + FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min) + FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) + CASE default ! Louis 1982 + zscf = SQRT(1.+CD*ABS(zri(i))) + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + END SELECT ! Calcul des drags - - - cdm(i)=cdmn(i)*FM(i) - cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) - - IF(nsrf.EQ.is_oce) THEN - - cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) - cdm(i)=MIN(cdm(i),cdmmax) - cdh(i)=MIN(cdh(i),cdhmax) - - ENDIF - - - ENDIF - -!xxxxxxxxxxx + cdm(i)=cdmn(i)*FM(i) + cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) + ENDIF + ENDIF ! fin du if (nsrf == is_oce) END DO ! Fin de la boucle sur l'horizontal -!xxxxxxxxxxx -! ================================================================= c END SUBROUTINE cdrag Index: LMDZ6/trunk/libf/phylmd/clc_core_cp.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/clc_core_cp.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/clc_core_cp.F90 (revision 4722) @@ -0,0 +1,232 @@ +! +! $Id$ +! + SUBROUTINE clc_core_cp(du, dt, dq, t1, q1, & + zu, zt, zq, & + P,& + n_it, mixte, & + coeffs, rugosm, rugosh) + + IMPLICIT NONE + + + ! INTEGER :: i,j,k + REAL, INTENT(IN) :: du,dt,dq,t1,q1 + REAL, INTENT(IN) :: zu,zt,zq, P + INTEGER, INTENT(IN):: n_it + LOGICAL, INTENT(IN) :: mixte + REAL, DIMENSION(3), INTENT(OUT) :: coeffs + real, intent(out) :: rugosm + real, intent(out) :: rugosh + + +! REAL :: du,dt,t1,dq,q1 + REAL :: dt_u, dq_u + + + + REAL :: cd, ch, cd_rt, cd_n10, ce_n10, ce,& + u_N, X,& + rho,cpa,le + + REAL :: usr,tsr,qsr + REAL, DIMENSION(3) :: zeta ! dans l'ordre: it courant, it precedente, var + REAL, DIMENSION(2,3) :: psi + REAL, PARAMETER :: g = 9.81, k=.41, sqrt_k = .640312 !avec g= gravité, k = cste von karman, sqrt_k = racine(k) + REAL :: logzu_10,logzt_10,logzq_10,logzt_zu,logzq_zu,& + u,cd_n10_rt,ch_n10_rt,ch_n10,chi,z_0m,z_0h,tv + INTEGER :: i,j,m1,m2,m3,m4,m5 + + + REAL :: phih,phid,phie + + + REAL, parameter :: alpha=2.7e-3,& !alpha = a1 ds formulation smith + beta = 1.42e-4,& !beta = a2 ds la formulation smith + gamma = 7.64e-5,& !Large and yaeger 2006 !gamma = a3 ds formulation smith +! gamma = 13.09e-3,& !Large and yaeger 2004 + q0 = 1.64474 ! q0 dzfinis mais pas reutilisée ds la suite ?? + + REAL, dimension(3) :: z + integer, dimension(2) :: tmp_shape + z = [zu, zt, zq] + + + + logzu_10 = LOG(zu/10.) + logzt_10 = LOG(zt/10.) + logzq_10 = LOG(zq/10.) + logzt_zu = logzt_10 - logzu_10 + logzq_zu = logzq_10 - logzu_10 + + cpa = 1004.67 !!! pas reutilise ds nvelle formulation... + + tv = t1 * (1+.608*q1) !!!!! tv pas reutiliser ds nvelle formulation... + rho = p / (287.1 * t1 * (1.+.61*q1)) !!!! rho pas reutiliser ds nvelle formulation.... + + + le = (2.501-.00237*(t1-273.15-dt))*1.e6 !!!!!! le pas reutiliser ds nvelle formulation.... + + u = max(du,.5) !!!! u = vitesse du vent 1 er niveau + u_N = u !!!!!!! u_N utilisée pour le calcul des premier calcul... + + + !!!!! initialisation des parametres et premier calcul de z_0 cd et ch + + + if (mixte) then + z_0m = 1.e-4 !!!!! pour l'initialisation on fixe z_0 a 10-4 + cd_n10 = k**2 / (log(10/z_0m))**2 !!!!! on calcul ensuite cd a 10 m d'après smith + cd_n10_rt = sqrt(cd_n10) + else +! Large and yaeger 2006 + cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur +! Large and yaeger 2004 +! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) + cd_n10_rt = sqrt(cd_n10) !!!!!! + z_0m = 10. * exp(-k/cd_n10_rt) !!!!!!! on calcul le premier z_0 + endif + + + if ( dt .gt. 0. ) then + ch_n10 = 0.018 * cd_n10_rt !!!!!!! si regime stable calcul du ch +! z_0h = 10. * exp(-k/ch_n10) + else + ch_n10 = 0.0327 * cd_n10_rt !!!!!!! si regime instable calcul du ch +! z_0h = 10. * exp(-k/ch_n10) + endif + + ce_n10 = 3.46e-2 *cd_n10_rt !!!!!! calcul ensuite de ce + + + dt_u = dt !!! def + dq_u = dq !!! def + + cd = cd_n10 !!! def + ch = ch_n10 !!! def + ce = ce_n10 !!! def + + cd_rt = sqrt(cd) !!! def + + ! open(7,file="err_rel.dat",position="append") + + +!!!!!!!!! début des itérations.... + + do i = 1,n_it + + usr = cd_rt * u + tsr = ch / cd_rt * dt_u + qsr = ce / cd_rt * dq_u + + + zeta = k*g / (usr**2) * tsr / t1& + * z + + zeta = sign( min(abs(zeta),10.0), zeta ) + + +!!!! calcul du zeta pour connaitre le signe et donc en deduire le regime... + + + + + IF ( zeta(1) .GT. 0 ) THEN !!! si regime stable alors valeurs pr chi et psi + + chi = 0.018 + + + psi(1,1) = -5.*zeta(1) + psi(2,1) = psi(1,1) + + ELSE !!!!! si regime instable alors valeurs pr psi et chi !!! + + chi = 0.0327 + + X = sqrt( sqrt( 1-16.*zeta(1) )) + psi(1,1) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) + psi(2,1) = 2.*LOG( (1 + X**2) / 2 ) + + END IF + + do j =2,3 + IF ( zeta(j) .GT. 0 ) THEN + psi(1,j) = -5.*zeta(j) + psi(2,j) = psi(1,1) + + ELSE + X = sqrt( sqrt( 1-16.*zeta(j) )) + psi(1,j) = 3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X) + psi(2,j) = 2.*LOG( (1 + X**2) / 2 ) + + END IF + enddo + + + dt_u = dt - tsr / k * & + ( logzt_zu + psi(2,1) - psi(2,2) ) + + dq_u = dq - qsr / k * & + ( logzq_zu + psi(2,1) - psi(2,3) ) + + + + u_N = u / ( 1 + cd_n10_rt / k * ( logzu_10 - psi(1,1)) ) + + if (mixte) then + z_0m=0.018*usr**2/9.81 + 0.11*14e-6 / (usr) +! z_0h=(0.11*14e-6 / (usr)) + 1.4e-5 + cd_n10 = k**2 / (log(10/z_0m))**2 + cd_n10_rt = sqrt(cd_n10) +! ch_n10 = k**2 / ((log(10/z_0m))*(log(10/z_0h))) +! ch_n10_rt = sqrt(ch_n10) + ch_n10 = 1.e-3 + else +! Large and yaeger 2006 + cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur +! Large and yaeger 2004 +! cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) + ! cd_n10_dun = -alpha / u_n**2 + gamma + cd_n10_rt = sqrt(cd_n10) + z_0m = 10. * exp(-k/cd_n10_rt) !!!!! c'est sur cette ligne la qu'il y a un bug... ou sur la ligne du u_N +! ch_n10 = chi * cd_n10_rt +! z_0h = 10. * exp(-k/ch_n10) +! ch_n10_rt = sqrt(ch_n10) + ch_n10 = chi * cd_n10_rt + endif + + + ce_n10 = 3.46e-2*cd_n10_rt + + + + phid = 1+cd_n10_rt/k * (logzu_10 - psi(1,1)) + phih = 1+chi/k * (logzu_10 - psi(2,1)) + phie = 1+3.46e-2/k * (logzu_10 - psi(2,1)) + + + cd_rt = cd_n10_rt / abs( phid ) + ch = ch_n10 / ( phih * phid ) + ce = ce_n10 / ( phie * phid ) + + cd=cd_rt**2 + + END DO + + usr = cd_rt * u + tsr = ch / cd_rt * dt_u + qsr = ce / cd_rt * dq_u + + + ! coeffs(:,m) = [ & + ! rho * usr**2,& + ! rho * cpa * usr *tsr,& + ! rho * le * usr *qsr] + + coeffs = [ cd_rt**2, ch , ce ] + + rugosm = z_0m + rugosh = z_0h + + RETURN + END subroutine clc_core_cp Index: LMDZ6/trunk/libf/phylmd/clesphys.h =================================================================== --- LMDZ6/trunk/libf/phylmd/clesphys.h (revision 4720) +++ LMDZ6/trunk/libf/phylmd/clesphys.h (revision 4722) @@ -43,4 +43,6 @@ !IM seuils cdrm, cdrh REAL cdmmax, cdhmax +!IM pour les params différentes Olivier Torres + INTEGER choix_bulk, nit_bulk, kz0 !IM param. stabilite s/ terres et en dehors REAL ksta, ksta_ter, f_ri_cd_min @@ -111,4 +113,5 @@ COMMON/clesphys/ & ! REAL FIRST +! rajout choix_bulk et nit_bulk kz0 par Olivier Torres & co2_ppm, solaire & & , RCO2, RCH4, RN2O, RCFC11, RCFC12 & @@ -138,4 +141,5 @@ & , ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad & & , iflag_con, nbapp_cv, nbapp_wk & + & , choix_bulk, nit_bulk, kz0 & & , iflag_ener_conserv & & , ok_suntime_rrtm & Index: LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90 (revision 4722) @@ -0,0 +1,631 @@ +!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SFX_LIC for details. version 1. +! ######### + +!------------------------ Modif Olivier Torres pour rajout fonction et pas appel ---------------- + + + +real function PSIFCTT(zet) + + real, intent(in) :: zet + + if(zet<0) then + x=(1.-(15*zet))**.5 + psik=2*log((1+x)/2) + x=(1.-(34.15*zet))**.3333 + psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) + f=zet*zet/(1+zet*zet) + PSIFCTT=(1-f)*psik+f*psic + + else + c=min(50.,.35*zet) + PSIFCTT=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) + endif +end FUNCTION PSIFCTT + +real function PSIFCTU(zet) + + real, intent(in) :: zet + + if (zet<0) then + + x=(1.-15.*zet)**.25 + psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) + x=(1.-10.15*zet)**.3333 + psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) + f=zet*zet/(1+zet*zet) + PSIFCTU=(1-f)*psik+f*psic + else + c=min(50.,.35*zet) + PSIFCTU=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) + endif +END FUNCTION PSIFCTU + +!----------------------------------- Fin Modif --------------------------------------------------- + + +SUBROUTINE COARE30_FLUX_CNRM(PZ0SEA,PTA,PSST,PQA, & + PVMOD,PZREF,PUREF,PPS,PQSATA,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,& + PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs ) +! ####################################################################### +! +! +!!**** *COARE25_FLUX* +!! +!! PURPOSE +!! ------- +! Calculate the surface fluxes of heat, moisture, and momentum over +! sea surface with bulk algorithm COARE3.0. +! +!!** METHOD +!! ------ +! transfer coefficients were obtained using a dataset which combined COARE +! data with those from three other ETL field experiments, and reanalysis of +! the HEXMAX data (DeCosmos et al. 1996). +! ITERMAX=3 +! Take account of the surface gravity waves on the velocity roughness and +! hence the momentum transfer coefficient +! NGRVWAVES=0 no gravity waves action (Charnock) !default value +! NGRVWAVES=1 wave age parameterization of Oost et al. 2002 +! NGRVWAVES=2 model of Taylor and Yelland 2001 +! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! Fairall et al (2003), J. of Climate, vol. 16, 571-591 +!! Fairall et al (1996), JGR, 3747-3764 +!! Gosnell et al (1995), JGR, 437-442 +!! Fairall et al (1996), JGR, 1295-1308 +!! +!! AUTHOR +!! ------ +!! C. Lebeaupin *Météo-France* (adapted from C. Fairall's code) +!! +!! MODIFICATIONS +!! ------------- +!! Original 1/06/2006 +!! B. Decharme 06/2009 limitation of Ri +!! B. Decharme 09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90 +!! B. Decharme 06/2013 bug in z0 (output) computation +!! J.Escobar 06/2013 for REAL4/8 add EPSILON management +!! C. Lebeaupin 03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value +!! add abort if no convergence +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +! +!USE MODD_SEAFLUX_n, ONLY : SEAFLUX_t +! + +!----------Rajout Olive --------- +USE dimphy +USE indice_sol_mod + +!-------------------------------- + +USE MODD_CSTS, ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, & + XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, & + XP00 + +!USE MODD_SURF_ATM, ONLY : XVZ0CM +! +!USE MODD_SURF_PAR, ONLY : XUNDEF, XSURF_EPSILON +!USE MODD_WATER_PAR +! +!USE MODI_SURFACE_RI +!USE MODI_WIND_THRESHOLD +!USE MODE_COARE30_PSI +! +!USE MODE_THERMOS +! +! +!USE MODI_ABOR1_SFX +! +! +!USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +!USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +! +! +!TYPE(SEAFLUX_t), INTENT(INOUT) :: S +! +REAL, DIMENSION(klon), INTENT(IN) :: PTA ! air temperature at atm. level (K) +REAL, DIMENSION(klon), INTENT(IN) :: PQA ! air humidity at atm. level (kg/kg) +!REAL, DIMENSION(:), INTENT(IN) :: PEXNA ! Exner function at atm. level +!REAL, DIMENSION(:), INTENT(IN) :: PRHOA ! air density at atm. level +REAL, DIMENSION(klon), INTENT(IN) :: PVMOD ! module of wind at atm. wind level (m/s) +REAL, DIMENSION(klon), INTENT(IN) :: PZREF ! atm. level for temp. and humidity (m) +REAL, DIMENSION(klon), INTENT(IN) :: PUREF ! atm. level for wind (m) +REAL, DIMENSION(klon), INTENT(IN) :: PSST ! Sea Surface Temperature (K) +!REAL, DIMENSION(:), INTENT(IN) :: PEXNS ! Exner function at sea surface +REAL, DIMENSION(klon), INTENT(IN) :: PPS ! air pressure at sea surface (Pa) +REAL, DIMENSION(klon), INTENT(IN) :: PRAIN !precipitation rate (kg/s/m2) +REAL, DIMENSION(klon), INTENT(IN) :: PPA ! air pressure at atm level (Pa) +REAL, DIMENSION(klon), INTENT(IN) :: PQSATA ! air pressure at atm level (Pa) +! +REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA! roughness length over the ocean +! +! surface fluxes : latent heat, sensible heat, friction fluxes +REAL, DIMENSION(klon), INTENT(OUT) :: PSFTH ! heat flux (W/m2) +REAL, DIMENSION(klon), INTENT(OUT) :: PSFTQ ! water flux (kg/m2/s) +REAL, DIMENSION(klon), INTENT(OUT) :: PUSTAR! friction velocity (m/s) +! +! diagnostics +REAL, DIMENSION(klon), INTENT(OUT) :: PQSAT ! humidity at saturation +REAL, DIMENSION(klon), INTENT(OUT) :: PCD ! heat drag coefficient +REAL, DIMENSION(klon), INTENT(OUT) :: PCDN ! momentum drag coefficient +REAL, DIMENSION(klon), INTENT(OUT) :: PCH ! neutral momentum drag coefficient +REAL, DIMENSION(klon), INTENT(OUT) :: PCE !transfer coef. for latent heat flux +REAL, DIMENSION(klon), INTENT(OUT) :: PRI ! Richardson number +REAL, DIMENSION(klon), INTENT(OUT) :: PRESA ! aerodynamical resistance +REAL, DIMENSION(klon), INTENT(OUT) :: PZ0HSEA ! heat roughness length + +LOGICAL, INTENT(IN) :: LPRECIP ! +LOGICAL, INTENT(IN) :: LPWG ! +real, dimension(3), intent(inout) :: coeffs +! +! + +!INCLUDE "YOMCST.h" +!INCLUDE "clesphys.h" + +!* 0.2 declarations of local variables +! +REAL, DIMENSION(SIZE(PTA)) :: ZVMOD ! wind intensity +REAL, DIMENSION(SIZE(PTA)) :: ZPA ! Pressure at atm. level +REAL, DIMENSION(SIZE(PTA)) :: ZTA ! Temperature at atm. level +REAL, DIMENSION(SIZE(PTA)) :: ZQASAT ! specific humidity at saturation at atm. level (kg/kg) +! +!rajout +REAL, DIMENSION(SIZE(PTA)) :: PEXNA ! Exner function at atm level +REAL, DIMENSION(SIZE(PTA)) :: PEXNS ! Exner function at atm level +! +REAL, DIMENSION(SIZE(PTA)) :: ZO ! rougness length ref +REAL, DIMENSION(SIZE(PTA)) :: ZWG ! gustiness factor (m/s) +! +REAL, DIMENSION(SIZE(PTA)) :: ZDU,ZDT,ZDQ,ZDUWG !differences +! +REAL, DIMENSION(SIZE(PTA)) :: ZUSR !velocity scaling parameter "ustar" (m/s) = friction velocity +REAL, DIMENSION(SIZE(PTA)) :: ZTSR !temperature sacling parameter "tstar" (degC) +REAL, DIMENSION(SIZE(PTA)) :: ZQSR !humidity scaling parameter "qstar" (kg/kg) +! +REAL, DIMENSION(SIZE(PTA)) :: ZU10,ZT10 !vertical profils (10-m height) +REAL, DIMENSION(SIZE(PTA)) :: ZVISA !kinematic viscosity of dry air +REAL, DIMENSION(SIZE(PTA)) :: ZO10,ZOT10 !roughness length at 10m +REAL, DIMENSION(SIZE(PTA)) :: ZCD,ZCT,ZCC +REAL, DIMENSION(SIZE(PTA)) :: ZCD10,ZCT10 !transfer coef. at 10m +REAL, DIMENSION(SIZE(PTA)) :: ZRIBU,ZRIBCU +REAL, DIMENSION(SIZE(PTA)) :: ZETU,ZL10 +! +REAL, DIMENSION(SIZE(PTA)) :: ZCHARN !Charnock number depends on wind module +REAL, DIMENSION(SIZE(PTA)) :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact +! +REAL, DIMENSION(SIZE(PTA)) :: ZZL,ZZTL!,ZZQL !Obukhovs stability + !param. z/l for u,T,q +REAL, DIMENSION(SIZE(PTA)) :: ZRR +REAL, DIMENSION(SIZE(PTA)) :: ZOT,ZOQ !rougness length ref +REAL, DIMENSION(SIZE(PTA)) :: ZPUZ,ZPTZ,ZPQZ !PHI funct. for u,T,q +! +REAL, DIMENSION(SIZE(PTA)) :: ZBF !constants to compute gustiness factor +! +REAL, DIMENSION(SIZE(PTA)) :: ZTAU !momentum flux (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZHF !sensible heat flux (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZEF !latent heat flux (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZWBAR !diag for webb correction but not used here after +REAL, DIMENSION(SIZE(PTA)) :: ZTAUR !momentum flux due to rain (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZRF !sensible heat flux due to rain (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZCHN,ZCEN !neutral coef. for heat and vapor +! +REAL, DIMENSION(SIZE(PTA)) :: ZLV !latent heat constant +! +REAL, DIMENSION(SIZE(PTA)) :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact +REAL, DIMENSION(SIZE(PTA)) :: ZXLR ! vaporisation heat at a given temperature +REAL, DIMENSION(SIZE(PTA)) :: ZCPLW ! specific heat for water at a given temperature +! +REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2 ! square of friction velocity +! +REAL, DIMENSION(SIZE(PTA)) :: ZDIRCOSZW! orography slope cosine (=1 on water!) +REAL, DIMENSION(SIZE(PTA)) :: ZAC ! Aerodynamical conductance +! +! +INTEGER, DIMENSION(SIZE(PTA)) :: ITERMAX ! maximum number of iterations +! +REAL :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants +REAL :: ZBETAGUST !gustiness factor +REAL :: ZZBL !atm. boundary layer depth (m) +REAL :: ZVISW !m2/s kinematic viscosity of water +REAL :: ZS !height of rougness length ref +REAL :: ZCH10 !transfer coef. at 10m + +REAL :: QSAT_SEAWATER +REAL :: QSATSEAW_1D +REAL, EXTERNAL :: PSIFCTU, PSIFCTT +! +INTEGER :: J, JLOOP !loop indice +!REAL(KIND=JPRB) :: ZHOOK_HANDLE + +!--------- Modif Olive ----------------- +REAL, DIMENSION(SIZE(PTA)) :: PRHOA +REAL, PARAMETER :: XUNDEF = 1.E+20 + +REAL :: XVCHRNK = 0.021 +REAL :: XVZ0CM = 1.0E-5 +!REAL :: XRIMAX + + +INTEGER :: PREF ! reference pressure for exner function +INTEGER :: NGRVWAVES ! Pour le choix du z0 + +INCLUDE "YOMCST.h" +INCLUDE "clesphys.h" + +!-------------------------------------- + + +PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:))) + +PREF = 100000. ! = 1000 hPa +NGRVWAVES = 1 + +PEXNA = (PPA/PREF)**(RD/RCPD) +PEXNS = (PPS/PREF)**(RD/RCPD) + +! +!------------------------------------------------------------------------------- +! +! 1. Initializations +! --------------- +! +! 1.1 Constants and parameters +! +!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',0,ZHOOK_HANDLE) +! +ZRVSRDM1 = XRV/XRD-1. ! 0.607766 +ZRDSRV = XRD/XRV ! 0.62198 +ZR2 = 1.-ZRDSRV ! pas utilisé dans cette routine +ZBETAGUST = 1.2 ! value based on TOGA-COARE experiment +ZZBL = 600. ! Set a default value for boundary layer depth +ZS = 10. ! Standard heigth =10m +ZCH10 = 0.00115 +! +ZVISW = 1.E-6 +! +! 1.2 Array initialization by undefined values +! +PSFTH (:)=XUNDEF +PSFTQ (:)=XUNDEF +PUSTAR(:)=XUNDEF +! +PCD(:) = XUNDEF +PCDN(:) = XUNDEF +PCH(:) = XUNDEF +PCE(:) =XUNDEF +PRI(:) = XUNDEF +! +PRESA(:)=XUNDEF +! +!------------------------------------------------------------------------------- +! 2. INITIAL GUESS FOR THE ITERATIVE METHOD +! ------------------------------------- +! +! 2.0 Temperature +! +! Set a non-zero value for the temperature gradient +! +WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.) + ZTA(:)=PTA(:)-1E-3 +ELSEWHERE + ZTA(:)=PTA(:) +ENDWHERE + +! 2.1 Wind and humidity +! +! Sea surface specific humidity +! +!PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:)) +PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:)) + + + +! +! Set a minimum value to wind +! +!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:)) + + +ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) ) !set a minimum value to wind +!ZVMOD = PVMOD !set a minimum value to wind + +! +! Specific humidity at saturation at the atm. level +! +ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD)) +!ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:)) +ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:)) + + +! +! +ZO(:) = 0.0001 +ZWG(:) = 0. +IF (LPWG) ZWG(:) = 0.5 +! +ZCHARN(:) = 0.011 +! +DO J=1,SIZE(PTA) + ! + ! 2.2 initial guess + ! + ZDU(J) = ZVMOD(J) !wind speed difference with surface current(=0) (m/s) + !initial guess for gustiness factor + ZDT(J) = -(ZTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference + ZDQ(J) = PQSAT(J)-PQA(J) !specific humidity difference + ! + ZDUWG(J) = SQRT(ZDU(J)**2+ZWG(J)**2) !wind speed difference including gustiness ZWG + ! + ! 2.3 initialization of neutral coefficients + ! + ZU10(J) = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J)) + ZUSR(J) = 0.035*ZU10(J) + ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+& + 8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11 + ! + ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J) + ZCD(J) = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2 !drag coefficient + ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2 + ZCT10(J)= ZCH10/SQRT(ZCD10(J)) + ZOT10(J)= ZS/EXP(XKARMAN/ZCT10(J)) + ! + !------------------------------------------------------------------------------- + ! Grachev and Fairall (JAM, 1997) + ZCT(J) = XKARMAN/LOG(PZREF(J)/ZOT10(J)) !temperature transfer coefficient + ZCC(J) = XKARMAN*ZCT(J)/ZCD(J) !z/L vs Rib linear coef. + ! + ZRIBCU(J) = -PUREF(J)/(ZZBL*0.004*ZBETAGUST**3) !saturation or plateau Rib + !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(ZTA(J)-XTT)*ZDQ)/& + ! &((ZTA(J)-XTT)*ZDUWG(J)**2) + ZRIBU(J) = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/& + (ZTA(J)*ZDUWG(J)**2) + ! + IF (ZRIBU(J)<0.) THEN + ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J)) !Unstable G and F + ELSE + ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable + ENDIF + ! + ZL10(J) = PUREF(J)/ZETU(J) !MO length + ! +ENDDO +! +! First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL) +ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF/ZL10)) +ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10)) +ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10)) +! +ZZL(:) = 0.0 +! +DO J=1,SIZE(PTA) + ! + IF (ZETU(J)>50.) THEN + ITERMAX(J) = 1 + ELSE + ITERMAX(J) = 3 !number of iterations + ENDIF + ! + !then modify Charnork for high wind speeds Chris Fairall's data + IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18.-10.) + IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018 + ! + ! 3. ITERATIVE LOOP TO COMPUTE USR, TSR, QSR + ! ------------------------------------------- + ! + ZHWAVE(J) = 0.018*ZVMOD(J)*ZVMOD(J)*(1.+0.015*ZVMOD(J)) + ZTWAVE(J) = 0.729*ZVMOD(J) + ZCWAVE(J) = XG*ZTWAVE(J)/(2.*XPI) + ZLWAVE(J) = ZTWAVE(J)*ZCWAVE(J) + ! +ENDDO +! + +! +DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop + ! + DO J=1,SIZE(PTA) + ! + IF (JLOOP.GT.ITERMAX(J)) CYCLE + ! + IF (NGRVWAVES==0) THEN + ZO(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG + 0.11*ZVISA(J)/ZUSR(J) !Smith 1988 + ELSE IF (NGRVWAVES==1) THEN + ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 & + + 0.11*ZVISA(J)/ZUSR(J) !Oost et al. 2002 + ELSE IF (NGRVWAVES==2) THEN + ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 & + + 0.11*ZVISA(J)/ZUSR(J) !Taulor and Yelland 2001 + ENDIF + ! + ZRR(J) = ZO(J)*ZUSR(J)/ZVISA(J) + ZOQ(J) = MIN(1.15E-4 , 5.5E-5/ZRR(J)**0.6) + ZOT(J) = ZOQ(J) + ! + ZZL(J) = XKARMAN * XG * PUREF(J) * & + ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / & + ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) ) + ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J) ! for T +! ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J) ! for Q + ENDDO + ! + ZPUZ(:) = PSIFCTU(ZZL(:)) + ZPTZ(:) = PSIFCTT(ZZTL(:)) + ! + DO J=1,SIZE(PTA) + ! + ! ZPQZ(J)=PSIFCTT(ZZQL(J)) + ZPQZ(J) = ZPTZ(J) + ! + ! 3.1 scale parameters + ! + ZUSR(J) = ZDUWG(J)*XKARMAN/(LOG(PUREF(J)/ZO(J)) -ZPUZ(J)) + ZTSR(J) = -ZDT(J) *XKARMAN/(LOG(PZREF(J)/ZOT(J))-ZPTZ(J)) + ZQSR(J) = -ZDQ(J) *XKARMAN/(LOG(PZREF(J)/ZOQ(J))-ZPQZ(J)) + ! + ! 3.2 Gustiness factor (ZWG) + ! + IF(LPWG) THEN + ZBF(J) = -XG/ZTA(J)*ZUSR(J)*(ZTSR(J)+ZRVSRDM1*ZTA(J)*ZQSR(J)) + IF (ZBF(J)>0.) THEN + ZWG(J) = ZBETAGUST*(ZBF(J)*ZZBL)**(1./3.) + ELSE + ZWG(J) = 0.2 + ENDIF + ENDIF + ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2) + ! + ENDDO + ! +ENDDO +!------------------------------------------------------------------------------- +! +! 4. COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES +! -------------------------------------------------------------- +! +ZTAU(:) = XUNDEF +ZHF(:) = XUNDEF +ZEF(:) = XUNDEF +! +ZWBAR(:) = 0. +ZTAUR(:) = 0. +ZRF(:) = 0. +! +DO J=1,SIZE(PTA) + ! + ! + ! 4. transfert coefficients PCD, PCH and PCE + ! and neutral PCDN, ZCHN, ZCEN + ! + PCD(J) = (ZUSR(J)/ZDUWG(J))**2. + PCH(J) = ZUSR(J)*ZTSR(J)/(ZDUWG(J)*(ZTA(J)*PEXNS(J)/PEXNA(J)-PSST(J))) + PCE(J) = ZUSR(J)*ZQSR(J)/(ZDUWG(J)*(PQA(J)-PQSAT(J))) + ! + PCDN(J) = (XKARMAN/LOG(ZS/ZO(J)))**2. + ZCHN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOT(J))) + ZCEN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOQ(J))) + ! + ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT) + ! + ! 4. 2 surface fluxes + ! + IF (ABS(PCDN(J))>1.E-2) THEN !!!! secure COARE3.0 CODE + write(*,*) 'pb PCDN in COARE30: ',PCDN(J) + write(*,*) 'point: ',J,"/",SIZE(PTA) + write(*,*) 'roughness: ', ZO(J) + write(*,*) 'ustar: ',ZUSR(J) + write(*,*) 'wind: ',ZDUWG(J) + CALL abort_physic('COARE30',': PCDN too large -> no convergence',1) + ELSE + ZTSR(J) = -ZTSR(J) + ZQSR(J) = -ZQSR(J) + ZTAU(J) = -PRHOA(J)*ZUSR(J)*ZUSR(J)*ZVMOD(J)/ZDUWG(J) + ZHF(J) = PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J) + ZEF(J) = PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J) + ! + ! 4.3 Contributions to surface fluxes due to rainfall + ! + ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour + ! adapter la formule de Clausius-Clapeyron (pour l'air + ! sec) au cas humide. + IF (LPRECIP) THEN + ! + ! heat surface fluxes + ! + ZTAC(J) = ZTA(J)-XTT + ! + ZXLR(J) = XLVTT + (XCPV-XCL)* ZTAC(J) ! latent heat of rain vaporization + ZDQSDT(J)= ZQASAT(J) * ZXLR(J) / (XRD*ZTA(J)**2) ! Clausius-Clapeyron relation + ZDTMP(J) = (1.0 + 3.309e-3*ZTAC(J) -1.44e-6*ZTAC(J)*ZTAC(J)) * & !heat diffusivity + 0.02411 / (PRHOA(J)*XCPD) + ! + ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94 ! water vapour diffusivity from eq (13.3) + ! ! of Pruppacher and Klett (1978) + ZALFAC(J)= 1.0 / (1.0 + & ! Eq.11 in GoF95 + ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD)) ! ZALFAC=wet-bulb factor (sans dim) + ZCPLW(J) = 4224.8482 + ZTAC(J) * & + ( -4.707 + ZTAC(J) * & + (0.08499 + ZTAC(J) * & + (1.2826e-3 + ZTAC(J) * & + (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat + ! + ZRF(J) = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * & !Eq.12 in GoF95 !SIGNE? + (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD ) + ! + ! Momentum flux due to rainfall + ! + ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96 + ! + ENDIF + ! + ! 4.4 Webb correction to latent heat flux + ! + ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) & + - ZUSR(J)*ZTSR(J)/ZTA(J) ! Eq.21*rhoa in FBR96 + ! + ! 4.5 friction velocity which contains correction du to rain + ! + ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J) + PUSTAR(J) = SQRT(ZUSTAR2(J)) + ! + ! 4.6 Total surface fluxes + ! + PSFTH (J) = ZHF(J) + ZRF(J) + PSFTQ (J) = ZEF(J) / ZLV(J) + ! + ENDIF +ENDDO + + +coeffs = [PCD,& + PCE,& + PCH] + +!------------------------------------------------------------------------------- +! +! 5. FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS +! ----------- +! 5.1 Richardson number +! +! +!------------STOP LA -------------------- +!ZDIRCOSZW(:) = 1. +! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,& +! PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI ) +!! +!! 5.2 Aerodynamical conductance and resistance +!! +!ZAC(:) = PCH(:)*ZVMOD(:) +!PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON) +! +!! 5.3 Z0 and Z0H over sea +!! +!PZ0SEA(:) = ZCHARN(:) * ZUSTAR2(:) / XG + XVZ0CM * PCD(:) / PCDN(:) +!! +!!PZ0HSEA(:) = PZ0SEA(:) +!! +!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE COARE30_FLUX_CNRM Index: LMDZ6/trunk/libf/phylmd/coare_cp.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/coare_cp.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/coare_cp.F90 (revision 4722) @@ -0,0 +1,252 @@ +real function psit_30(zet) + + real, intent(in) :: zet + + if(zet<0) then + x=(1.-(15*zet))**.5 + psik=2*log((1+x)/2) + x=(1.-(34.15*zet))**.3333 + psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) + f=zet*zet/(1+zet*zet) + psit_30=(1-f)*psik+f*psic + + else + c=min(50.,.35*zet) + psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) + endif +end FUNCTION psit_30 + +real function psiuo(zet) + + real, intent(in) :: zet + + if (zet<0) then + + x=(1.-15.*zet)**.25 + psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) + x=(1.-10.15*zet)**.3333 + psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) + f=zet*zet/(1+zet*zet) + psiuo=(1-f)*psik+f*psic + else + c=min(50.,.35*zet) + psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) + endif +END FUNCTION psiuo + + +real function gen_q1(ts,p,dq) + + real, intent(in) :: ts,p,dq + real es,qs + + es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p) + qs = es*62.197/(p-37.8*es) + + q1 = dq + qs + +end function gen_q1 + + + + +subroutine coare_cp(du,dt,dq,& + t,q,& + zu,zt,zq,& + p,zi,& + nits,& + coeffs,rugosm,rugosh) + ! zo,tau,hsb,hlb,& + ! var) + !version with shortened iteration modified Rt and Rq + !uses wave information wave period in s and wave ht in m + !no wave, standard coare 2.6 charnock: jwave=0 + !Oost et al. zo=50/2/pi L (u*/c)**4.5 if jwave=1 + !taylor and yelland zo=1200 h*(L/h)**4.5 jwave=2 + ! + USE MODD_CSTS, ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, & + XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, & + XP00 + + IMPLICIT NONE + + real, intent(in) :: du,dt,dq,t,q + real, intent(in) :: zu,zt,zq,p,zi + integer, intent(in) :: nits + ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb + ! real, dimension(nits,3), intent(out) :: var + real, dimension(3), intent(out) :: coeffs + real, intent(out) :: rugosm + real, intent(out) :: rugosh + + real, parameter :: beta=1.2, von=.4, fdg = 1. ,& + tdk = 273.16, pi = 3.141593, grav = 9.82,& + rgas = 287.1 + + integer, dimension(3) :: shape_input + + real bf, cc, rhoa, visa,& + u10, ut, uts, ut0, ug,& + cd10, ch10, ct, ct10,& + charn, ribu,& + l, l10, & + ribcu,rr,& + zet,zetu, zom10, zoh10, zot, zoq, zot10,& + cd, ch, le, cpa, cpv,& + usr,qsr,tsr,& + zom, t_c + +!, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE + + + + + real old_usr, old_tsr, old_qsr,tmp + + real, external :: psit_30, psiuo, grv + integer i,j,k +!---------------- Rajout pour prendre en compte différent Z0 --------------------------------! +! INTEGER :: NGRVWAVES ! Pour le choix du z0 +! NGRVWAVES = 2 +!--------------------------------------------------------------------------------------------- +!-------------------- Attention Modif réalisée pas SURRR ------------------------------------- +!--------------------------------------------------------------------------------------------- + + Ribcu=-zu/zi/.004/Beta**3 + cpa=1004.67 + + t_c = t - 273.15 + + + Le=(2.501-.00237*(t_c-dt))*1e6 + + + visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) + + + cpv=cpa*(1+0.84*Q) + rhoa=P/(Rgas*t*(1+0.61*Q)) + + + + + + + ug=.2 + + ut=sqrt(du*du+ug*ug) + ut=MAX(ut , 0.1 * MIN(10.,zu) ) + + u10=ut*log(10/1e-4)/log(zu/1e-4) + usr=.035*u10 ! turbulent friction velocity (m/s), including gustiness + zom10=0.011*usr*usr/grav+0.11*visa/usr ! roughness length for u (smith 88) + Cd10=(von/log(10/zom10))**2 +! zoh10=0.40*visa/usr+1.4e-5 +! Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca + Ch10=0.00115 + Ct10=Ch10/sqrt(Cd10) + zot10=10/exp(von/Ct10) ! roughness length for t + Cd=(von/log(zu/zom10))**2 + Ct=von/log(zt/zot10) + CC=von*Ct/Cd + + ut0 = ut + + + ut = ut0 + Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 + + if (Ribu .LT. 0) then + zetu=CC*Ribu/(1+Ribu/Ribcu) + else + zetu=CC*Ribu*(1+27/9*Ribu/CC) + endif + + L10=zu/zetu + + ! if (zetu .GT. 50) then + ! nits=1 + ! endif + + usr=ut*von/(log(zu/zom10)-psiuo(zetu)) + tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) + qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) + ! tkt=.001 + + ! charnock constant - lin par morceau - constant + if ( ut <= 10. ) then + charn=0.011 + else + if (ut .GT. 18) then + charn=0.018 + else + charn=0.011+(ut-10)/(18-10)*(0.018-0.011) + endif + endif + +! ZHWAVE = 0.018*ut*ut*(1.+0.015*ut) +! ZTWAVE = 0.729*ut +! ZCWAVE = XG*ZTWAVE/(2.*pi) +! ZLWAVE = ZTWAVE*ZCWAVE + + + !*************** bulk loop ************ + do i=1, nits + + zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q) + +! IF (NGRVWAVES==0) THEN +! zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988 +! ELSE IF (NGRVWAVES==1) THEN +! zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 & +! + 0.11*visa/usr !Oost et al. 2002 +! ELSE IF (NGRVWAVES==2) THEN +! zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 & +! + 0.11*visa/usr !Taulor and Yelland 2001 +! ENDIF + + zom=charn*usr*usr/grav+0.11*visa/usr + + rr=zom*usr/visa + L=zu/zet + zoq=min(1.15e-4,5.5e-5/rr**.6) ! a modifier + zot=zoq ! a modifier +! zot=0.40*visa/usr+1.4e-5 + + old_usr = usr + old_tsr = tsr + old_qsr = qsr + + usr=ut*von/(log(zu/zom)-psiuo(zu/L)) + tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) + qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) + + Bf=-grav/t*usr*(tsr+.61*t*qsr) + + if (Bf .GT. 0) then + ug=Beta*(Bf*zi)**.333 + else + ug=.2 + endif + + ut=sqrt(du*du+ug*ug) + ut=MAX(ut , 0.1 * MIN(10.,zu) ) + + + enddo !bulk iter loop + + ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress + ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr + ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr + + tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du + + coeffs = [tmp**2,& + von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,& + von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp] + + rugosm = zom + rugosh = zot + + +end subroutine coare_cp Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 4722) @@ -243,4 +243,8 @@ LOGICAL, SAVE :: ok_new_lscp_omp LOGICAL, SAVE :: ok_icefra_lscp_omp + !rajout de choix_bulk et nit_bulk par Olivier Torres + INTEGER,SAVE :: choix_bulk_omp + INTEGER,SAVE :: nit_bulk_omp + INTEGER,SAVE :: kz0_omp LOGICAL, SAVE :: ok_bs_omp, ok_rad_bs_omp @@ -936,4 +940,31 @@ nbapp_rad_omp = 12 CALL getin('nbapp_rad',nbapp_rad_omp) + + !rajout Olivier Torres + !Config Key = choix_bulk + !Config Desc = choix de la formulation bulk a prendre dans clcdrag au-dessus de l'ocean + !Config Def = 0 + !Config 0 -> originale (lmdz/Louis 79) + !Config 1 -> COARE + !Config 2 -> CORE-"pure" (cf. Large) + !Config 3 -> CORE-"mixte" (avec z_0 et C_T^N donnees par Smith 88) + choix_bulk_omp = 0 + CALL getin('choix_bulk',choix_bulk_omp) + + !Config Key = nit_bulk + !Config Desc = choix du nombre d'it de pt fixe dans la bulk + !Config Def = 5 + nit_bulk_omp = 1 + CALL getin('nit_bulk',nit_bulk_omp) + + !Config Key = kz0 + !Config Desc = choix de la formulation z0 pour la bulk ECUME + !Config Def = 1 + !Config 0 -> ARPEGE formulation + !Config 1 -> Smith Formulation + !Config 2 -> Direct computation using the stability functions + kz0_omp = 0 + CALL getin('kz0',kz0_omp) + !Config Key = iflag_con @@ -2496,4 +2527,8 @@ var_fco2_land_cor = var_fco2_land_cor_omp dms_cycle_cpl = dms_cycle_cpl_omp + !rajout Olivier Torres + kz0=kz0_omp + choix_bulk = choix_bulk_omp + nit_bulk = nit_bulk_omp ! Test of coherence between type_ocean and version_ocean @@ -2847,4 +2882,8 @@ WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol + !rajout Olivier Torres + write(lunout,*) 'choix_bulk = ', choix_bulk + write(lunout,*) 'nit_bulk = ', nit_bulk + write(lunout,*) 'kz0 = ', kz0 !$OMP END MASTER Index: LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90 (revision 4722) @@ -0,0 +1,777 @@ +!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SFX_LIC for details. version 1. +! ######### + SUBROUTINE ECUMEV6_FLUX(PZ0SEA,PTA,PSST,PQA,PQSAT,PVMOD, & + PZREF,PSSS,PUREF,PPS,PPA,OPRECIP,OPWEBB, & + PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE, & + PRI,PRESA,PRAIN,PZ0HSEA,OPERTFLUX,coeffs ) +!############################################################################### +!! +!!**** *ECUMEV6_FLUX* +!! +!! PURPOSE +!! ------- +! Calculate the surface turbulent fluxes of heat, moisture, and momentum +! over sea surface + corrections due to rainfall & Webb effect. +!! +!!** METHOD +!! ------ +! The estimation of the transfer coefficients relies on the iterative +! computation of the scaling parameters U*/Teta*/q*. The convergence is +! supposed to be reached in NITERFL iterations maximum. +! Neutral transfer coefficients for momentum/temperature/humidity +! are computed as a function of the 10m-height neutral wind speed using +! the ECUME_V6 formulation based on the multi-campaign (POMME,FETCH,CATCH, +! SEMAPHORE,EQUALANT) ALBATROS dataset. +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! Fairall et al (1996), JGR, 3747-3764 +!! Gosnell et al (1995), JGR, 437-442 +!! Fairall et al (1996), JGR, 1295-1308 +!! +!! AUTHOR +!! ------ +!! C. Lebeaupin *Météo-France* (adapted from S. Belamari's code) +!! +!! MODIFICATIONS +!! ------------- +!! Original 15/03/2005 +!! Modified 01/2006 C. Lebeaupin (adapted from A. Pirani's code) +!! Modified 08/2009 B. Decharme: limitation of Ri +!! Modified 09/2012 B. Decharme: CD correction +!! Modified 09/2012 B. Decharme: limitation of Ri in surface_ri.F90 +!! Modified 10/2012 P. Le Moigne: extra inputs for FLake use +!! Modified 06/2013 B. Decharme: bug in z0 (output) computation +!! Modified 12/2013 S. Belamari: ZRF computation updated: +!! 1. ZP00/PPA in ZDWAT, ZLVA in ZDQSDT/ZBULB/ZRF +!! 2. ZDWAT/ZDTMP in ZBULB/ZRF (Gosnell et al 95) +!! 3. cool skin correction included +!! Modified 01/2014 S. Belamari: salinity impact on latent heat of +!! vaporization of seawater included +!! Modified 01/2014 S. Belamari: new formulation for pure water +!! specific heat (ZCPWA) +!! Modified 01/2014 S. Belamari: 4 choices for PZ0SEA computation +!! Modified 12/2015 S. Belamari: ECUME now provides parameterisations +!! for: U10n*sqrt(CDN) instead of CDN +!! U10n*CHN/sqrt(CDN) " CHN +!! U10n*CEN/sqrt(CDN) " CEN +!! Modified 01/2016 S. Belamari: New ECUME formulation +!! +!! To be done: +!! include gustiness computation following Mondon & Redelsperger (1998) +!!! +!------------------------------------------------------------------------------- +!! +!! MODIFICATIONS RELATED TO SST CORRECTION COMPUTATION +!! --------------------------------------------------- +!! Modified 09/2013 S. Belamari: use 0.98 for the ocean emissivity +!! following up to date satellite measurements in +!! the 8-14 μm range (obtained values range from +!! 0.98 to 0.99). +!!! +!------------------------------------------------------------------------------- +! +! 0. DECLARATIONS +! ------------ +! +USE dimphy +USE indice_sol_mod +USE MODD_CSTS, ONLY : XPI, XDAY, XKARMAN, XG, XP00, XSTEFAN, XRD, XRV, & + XCPD, XCPV, XCL, XTT, XLVTT + + +!USE MODD_SURF_PAR, ONLY : XUNDEF +!USE MODD_SURF_ATM, ONLY : XVCHRNK, XVZ0CM +!USE MODD_REPROD_OPER, ONLY : CCHARNOCK +! +!USE MODE_THERMOS +!USE MODI_WIND_THRESHOLD +!USE MODI_SURFACE_RI +! +!USE YOMHOOK, ONLY : LHOOK, DR_HOOK +!USE PARKIND1, ONLY : JPRB +! +!USE MODI_ABOR1_SFX +! +IMPLICIT NONE +! +! 0.1. Declarations of arguments +! +REAL, DIMENSION(klon), INTENT(IN) :: PVMOD ! module of wind at atm level (m/s) +REAL, DIMENSION(klon), INTENT(IN) :: PTA ! air temperature at atm level (K) +REAL, DIMENSION(klon), INTENT(IN) :: PQA ! air spec. hum. at atm level (kg/kg) +REAL, DIMENSION(klon), INTENT(IN) :: PQSAT ! sea surface spec. hum. (kg/kg) +REAL, DIMENSION(klon), INTENT(IN) :: PPA ! air pressure at atm level (Pa) +!REAL, DIMENSION(:), INTENT(IN) :: PRHOA ! air density at atm level (kg/m3) +!REAL, DIMENSION(:), INTENT(IN) :: PEXNA ! Exner function at atm level +REAL, DIMENSION(klon), INTENT(IN) :: PUREF ! atm level for wind (m) +REAL, DIMENSION(klon), INTENT(IN) :: PZREF ! atm level for temp./hum. (m) +REAL, DIMENSION(klon), INTENT(IN) :: PSSS ! Sea Surface Salinity (g/kg) +REAL, DIMENSION(klon), INTENT(IN) :: PPS ! air pressure at sea surface (Pa) +!REAL, DIMENSION(:), INTENT(IN) :: PEXNS ! Exner function at sea surface +!REAL, DIMENSION(:), INTENT(IN) :: PPERTFLUX ! stochastic flux perturbation pattern +! for correction +!REAL, INTENT(IN) :: PICHCE ! +LOGICAL, INTENT(IN) :: OPRECIP ! +LOGICAL, INTENT(IN) :: OPWEBB ! +LOGICAL, INTENT(IN) :: OPERTFLUX +REAL, DIMENSION(klon), INTENT(IN) :: PRAIN ! precipitation rate (kg/s/m2) +! +!INTEGER, INTENT(IN) :: KZ0 +! +REAL, DIMENSION(klon), INTENT(INOUT) :: PSST ! Sea Surface Temperature (K) +REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA ! roughness length over sea +REAL, DIMENSION(klon), INTENT(OUT) :: PZ0HSEA ! heat roughness length over sea + +! surface fluxes : latent heat, sensible heat, friction fluxes +REAL, DIMENSION(klon), INTENT(OUT) :: PUSTAR ! friction velocity (m/s) +REAL, DIMENSION(klon), INTENT(OUT) :: PSFTH ! heat flux (W/m2) +REAL, DIMENSION(klon), INTENT(OUT) :: PSFTQ ! water flux (kg/m2/s) + +! diagnostics +REAL, DIMENSION(klon), INTENT(OUT) :: PCD ! transfer coef. for momentum +REAL, DIMENSION(klon), INTENT(OUT) :: PCH ! transfer coef. for temperature +REAL, DIMENSION(klon), INTENT(OUT) :: PCE ! transfer coef. for humidity +REAL, DIMENSION(klon), INTENT(OUT) :: PCDN ! neutral coef. for momentum +REAL, DIMENSION(klon), INTENT(OUT) :: PRI ! Richardson number +REAL, DIMENSION(klon), INTENT(OUT) :: PRESA ! aerodynamical resistance +real, dimension(3), intent(out) :: coeffs + +! 0.2. Declarations of local variables +! +! specif SB +INTEGER, DIMENSION(SIZE(PTA)) :: JCV ! convergence index +INTEGER, DIMENSION(SIZE(PTA)) :: JITER ! nb of iterations to converge +!rajout +REAL, DIMENSION(SIZE(PTA)) :: PEXNA ! Exner function at atm level +REAL, DIMENSION(SIZE(PTA)) :: PEXNS ! Exner function at atm level +! +REAL, DIMENSION(SIZE(PTA)) :: ZTAU ! momentum flux (N/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZHF ! sensible heat flux (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZEF ! latent heat flux (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZTAUR ! momentum flx due to rain (N/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZRF ! sensible flx due to rain (W/m2) +REAL, DIMENSION(SIZE(PTA)) :: ZEFWEBB ! Webb corr. on latent flx (W/m2) + +REAL, DIMENSION(SIZE(PTA)) :: ZVMOD ! wind intensity at atm level (m/s) +REAL, DIMENSION(SIZE(PTA)) :: ZQSATA ! sat.spec.hum. at atm level (kg/kg) +REAL, DIMENSION(SIZE(PTA)) :: ZLVA ! vap.heat of pure water at atm level (J/kg) +REAL, DIMENSION(SIZE(PTA)) :: ZLVS ! vap.heat of seawater at sea surface (J/kg) +REAL, DIMENSION(SIZE(PTA)) :: ZCPA ! specif.heat moist air (J/kg/K) +REAL, DIMENSION(SIZE(PTA)) :: ZVISA ! kinemat.visc. of dry air (m2/s) +REAL, DIMENSION(SIZE(PTA)) :: ZDU ! U vert.grad. (real atm) +REAL, DIMENSION(SIZE(PTA)) :: ZDT,ZDQ ! T,Q vert.grad. (real atm) +REAL, DIMENSION(SIZE(PTA)) :: ZDDU ! U vert.grad. (real atm + gust) +REAL, DIMENSION(SIZE(PTA)) :: ZDDT,ZDDQ ! T,Q vert.grad. (real atm + WL/CS) +REAL, DIMENSION(SIZE(PTA)) :: ZUSR ! velocity scaling param. (m/s) + ! =friction velocity +REAL, DIMENSION(SIZE(PTA)) :: ZTSR ! temperature scaling param. (K) +REAL, DIMENSION(SIZE(PTA)) :: ZQSR ! humidity scaling param. (kg/kg) +REAL, DIMENSION(SIZE(PTA)) :: ZDELTAU10N,ZDELTAT10N,ZDELTAQ10N + ! U,T,Q vert.grad. (10m, neutral atm) +REAL, DIMENSION(SIZE(PTA)) :: ZUSR0,ZTSR0,ZQSR0 ! ITERATIVE PROCESS +REAL, DIMENSION(SIZE(PTA)) :: ZDUSTO,ZDTSTO,ZDQSTO ! ITERATIVE PROCESS +REAL, DIMENSION(SIZE(PTA)) :: ZPSIU,ZPSIT! PSI funct for U, T/Q (Z0 comp) +REAL, DIMENSION(SIZE(PTA)) :: ZCHARN ! Charnock parameter (Z0 comp) + +REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2 ! square of friction velocity +REAL, DIMENSION(SIZE(PTA)) :: ZAC ! aerodynamical conductance +REAL, DIMENSION(SIZE(PTA)) :: ZDIRCOSZW ! orography slope cosine + ! (=1 on water!) +REAL, DIMENSION(SIZE(PTA)) :: ZPARUN,ZPARTN,ZPARQN ! neutral parameter for U,T,Q + +!-- rajout pour la pression saturante +REAL, DIMENSION(SIZE(PPA)) :: ZFOES ! [OPWEBB] +REAL, DIMENSION(SIZE(PPA)) :: ZWORK1 +REAL, DIMENSION(SIZE(PPA)) :: ZWORK2 +REAL, DIMENSION(SIZE(PPS)) :: ZWORK1A +REAL, DIMENSION(SIZE(PPS)) :: ZWORK2A +!##################### + +REAL, DIMENSION(0:5) :: ZCOEFU,ZCOEFT,ZCOEFQ + +!--------- Modif Olive ----------------- +REAL, DIMENSION(SIZE(PTA)) :: PRHOA +REAL, PARAMETER :: XUNDEF = 1.E+20 + + +REAL :: XVCHRNK = 0.021 +REAL :: XVZ0CM = 1.0E-5 +!REAL :: XRIMAX + +CHARACTER :: CCHARNOCK = 'NEW' + + +!-------------------------------------- + + +! local constants +LOGICAL :: OPCVFLX ! to force convergence +INTEGER :: NITERMAX ! nb of iterations to get free convergence +INTEGER :: NITERSUP ! nb of additional iterations if OPCVFLX=.TRUE. +INTEGER :: NITERFL ! maximum number of iterations +REAL :: ZETV,ZRDSRV ! thermodynamic constants +REAL :: ZSQR3 +REAL :: ZLMOMIN,ZLMOMAX ! min/max value of Obukhovs stability param. z/l +REAL :: ZBTA,ZGMA ! parameters of the stability functions +REAL :: ZDUSR0,ZDTSR0,ZDQSR0 ! maximum gap for USR/TSR/QSR between 2 steps +REAL :: ZP00 ! [OPRECIP] - water vap. diffusiv.ref.press.(Pa) +REAL :: ZUTU,ZUTT,ZUTQ ! U10n threshold in ECUME parameterisation +REAL :: ZCDIRU,ZCDIRT,ZCDIRQ ! coef directeur pour fonction affine U,T,Q +REAL :: ZORDOU,ZORDOT,ZORDOQ ! ordonnee a l'origine pour fonction affine U,T,Q + +INTEGER :: JJ ! for ITERATIVE PROCESS +INTEGER :: JLON,JK +REAL :: ZLMOU,ZLMOT ! Obukhovs param. z/l for U, T/Q +REAL :: ZPSI_U,ZPSI_T ! PSI funct. for U, T/Q +REAL :: Z0TSEA,Z0QSEA ! roughness length for T, Q +REAL :: ZCHIC,ZCHIK,ZPSIC,ZPSIK,ZLOGUS10,ZLOGTS10 +REAL :: ZTAC,ZCPWA,ZDQSDT,ZDWAT,ZDTMP,ZBULB ! [OPRECIP] +REAL :: ZWW ! [OPWEBB] + + +INTEGER :: PREF ! reference pressure for exner function +REAL, DIMENSION(klon) :: PQSATA ! sea surface spec. hum. (kg/kg) + +REAL :: qsat_seawater2,qsat_seawater + +INCLUDE "YOMCST.h" +INCLUDE "clesphys.h" + +!REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------- +!----------------------- Modif Olive calcul de PRHOA --------------------------- + +!write(*,*) "PZ0SEA ",PZ0SEA +!write(*,*) "PTA ",PTA +!write(*,*) "PSST ",PSST +!write(*,*) "PQA ",PQA +!write(*,*) "PVMOD ",PVMOD +!write(*,*) "PZREF ",PZREF +!write(*,*) "PUREF ",PUREF +!write(*,*) "PPS ",PPS +!write(*,*) "PPA ",PPA +!write(*,*) "OPRECIP ",OPRECIP +!write(*,*) "PZ0HSEA ",PZ0HSEA +!write(*,*) "PRAIN ",PRAIN + + +PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:))) +!write(*,*) "klon klon ",klon,PTA +!write(*,*) "PRHOA ",SIZE(PRHOA),PRHOA + +PREF = 100900. ! = 1000 hPa + +!PEXNA = (PPA/PPS)**RKAPPA +!PEXNS = (PPS/PPS)**RKAPPA + +PEXNA = (PPA/PREF)**(RD/RCPD) +PEXNS = (PPS/PREF)**(RD/RCPD) + +!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',0,ZHOOK_HANDLE) +! +ZDUSR0 = 1.E-06 +ZDTSR0 = 1.E-06 +ZDQSR0 = 1.E-09 +! +NITERMAX = 5 +NITERSUP = 5 +OPCVFLX = .TRUE. +! +NITERFL = NITERMAX +IF (OPCVFLX) NITERFL = NITERMAX+NITERSUP +! +ZCOEFU = (/ 1.00E-03, 3.66E-02, -1.92E-03, 2.32E-04, -7.02E-06, 6.40E-08 /) +ZCOEFT = (/ 5.36E-03, 2.90E-02, -1.24E-03, 4.50E-04, -2.06E-05, 0.0 /) +ZCOEFQ = (/ 1.00E-03, 3.59E-02, -2.87E-04, 0.0, 0.0, 0.0 /) +! +ZUTU = 40.0 +ZUTT = 14.4 +ZUTQ = 10.0 +! +ZCDIRU = ZCOEFU(1) + 2.0*ZCOEFU(2)*ZUTU + 3.0*ZCOEFU(3)*ZUTU**2 & + + 4.0*ZCOEFU(4)*ZUTU**3 + 5.0*ZCOEFU(5)*ZUTU**4 +ZCDIRT = ZCOEFT(1) + 2.0*ZCOEFT(2)*ZUTT + 3.0*ZCOEFT(3)*ZUTT**2 & + + 4.0*ZCOEFT(4)*ZUTT**3 +ZCDIRQ = ZCOEFQ(1) + 2.0*ZCOEFQ(2)*ZUTQ +! +ZORDOU = ZCOEFU(0) + ZCOEFU(1)*ZUTU + ZCOEFU(2)*ZUTU**2 + ZCOEFU(3)*ZUTU**3 & + + ZCOEFU(4)*ZUTU**4 + ZCOEFU(5)*ZUTU**5 +ZORDOT = ZCOEFT(0) + ZCOEFT(1)*ZUTT + ZCOEFT(2)*ZUTT**2 + ZCOEFT(3)*ZUTT**3 & + + ZCOEFT(4)*ZUTT**4 +ZORDOQ = ZCOEFQ(0) + ZCOEFQ(1)*ZUTQ + ZCOEFQ(2)*ZUTQ**2 +! +!------------------------------------------------------------------------------- +! +! 1. AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES. +! -------------------------------------------------------------------- +! +ZDIRCOSZW(:) = 1.0 +! +ZETV = XRV/XRD-1.0 !~0.61 (cf Liu et al. 1979) +ZRDSRV = XRD/XRV !~0.622 +ZSQR3 = SQRT(3.0) +ZLMOMIN = -200.0 +ZLMOMAX = 0.25 +ZBTA = 16.0 +ZGMA = 7.0 !initially =4.7, modified to 7.0 following G. Caniaux +! +ZP00 = 1013.25E+02 +! +PCD = XUNDEF +PCH = XUNDEF +PCE = XUNDEF +PCDN = XUNDEF +ZUSR = XUNDEF +ZTSR = XUNDEF +ZQSR = XUNDEF +ZTAU = XUNDEF +ZHF = XUNDEF +ZEF = XUNDEF +! +PSFTH = XUNDEF +PSFTQ = XUNDEF +PUSTAR = XUNDEF +PRESA = XUNDEF +PRI = XUNDEF +! +ZTAUR = 0.0 +ZRF = 0.0 +ZEFWEBB = 0.0 +! +!------------------------------------------------------------------------------- +! +! 2. INITIALISATIONS BEFORE ITERATIVE LOOP. +! ------------------------------------------- +! +!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:)) !set a minimum value to wind +ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) ) !set a minimum value to wind + +write(*,*) "ZVMOD ",SIZE(ZVMOD) + +! +! 2.0. Radiative fluxes - For warm layer & cool skin +! +! 2.0b. Warm Layer correction +! +! 2.1. Specific humidity at saturation +! +WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF) + PQSATA = QSAT_SEAWATER2 (PSST(:),PPS(:),PSSS(:)) !at sea surface +ELSEWHERE + PQSATA (:) = QSAT_SEAWATER (PSST(:),PPS(:)) !at sea surface +ENDWHERE + +!ZQSATA(:) = QSAT(PTA(:),PPA(:)) !at atm level + +!### OLIVIER POUR PRESSION SATURANTE ##### +!------------------------------------------------------------------------------- +! +ZFOES = 1 !PSAT(PT(:)) +ZFOES = 0.98*ZFOES +! +ZWORK1 = ZFOES/PPS +ZWORK2 = XRD/XRV + +ZWORK1A = ZFOES/PPA +ZWORK2A = XRD/XRV + +!write(*,*) "ZFOES ",ZFOES +!write(*,*) "PPS ",PPS +!write(*,*) "ZWORK1 ",ZWORK1 +!write(*,*) "XRD ",XRD +!write(*,*) "XRV ",XRV +!write(*,*) "PPA ",PPA +!write(*,*) "ZWORK1A ",ZWORK1A + +write(*,*) "PQSAT : ",PQSAT +write(*,*) "PQSATA : ",PQSATA + + +! +!* 2. COMPUTE SATURATION HUMIDITY +! --------------------------- +! +!PQSAT = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1) +!ZQSATA = ZWORK2A*ZWORK1A / (1.+(ZWORK2A-1.)*ZWORK1A) +ZQSATA = QSAT_SEAWATER (PTA(:),PPA(:)) !at sea surface + + +! +! 2.2. Gradients at the air-sea interface +! +ZDU(:) = ZVMOD(:) !one assumes u is measured / sea surface current +ZDT(:) = PTA(:)/PEXNA(:)-PSST(:)/PEXNS(:) +ZDQ(:) = PQA(:)-PQSATA(:) + +write(*,*) "PQA ",PQA(:) +write(*,*) "PQSAT",PQSAT(:) +write(*,*) "ZDQ",ZDQ(:) +! +! 2.3. Latent heat of vaporisation +! +ZLVA(:) = XLVTT+(XCPV-XCL)*(PTA (:)-XTT) !of pure water at atm level +ZLVS(:) = XLVTT+(XCPV-XCL)*(PSST(:)-XTT) !of pure water at sea surface + + + +write(*,*) "ZLVA ",ZLVA +write(*,*) "ZLVS ",ZLVS + + +WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF) + ZLVS(:) = ZLVS(:)*(1.0-1.00472E-3*PSSS(:)) !of seawater at sea surface +ENDWHERE +! +! 2.4. Specific heat of moist air (Businger 1982) +! +!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:)) +ZCPA(:) = XCPD +! +! 2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11) +! +ZVISA(:) = 1.326E-05*(1.0+6.542E-03*(PTA(:)-XTT)+8.301E-06*(PTA(:)-XTT)**2 & + -4.84E-09*(PTA(:)-XTT)**3) +! +! 2.4c Coefficients for warm layer and/or cool skin correction +! +! 2.5. Initial guess +! +ZDDU(:) = ZDU(:) +ZDDT(:) = ZDT(:) +ZDDQ(:) = ZDQ(:) +ZDDU(:) = SIGN(MAX(ABS(ZDDU(:)),10.0*ZDUSR0),ZDDU(:)) +ZDDT(:) = SIGN(MAX(ABS(ZDDT(:)),10.0*ZDTSR0),ZDDT(:)) +ZDDQ(:) = SIGN(MAX(ABS(ZDDQ(:)),10.0*ZDQSR0),ZDDQ(:)) + +write(*,*) "ZDDU ",ZDDU +write(*,*) "ZDDQ ",ZDDQ +write(*,*) "ZDDT ",ZDDT + +! +JCV (:) = -1 +ZUSR(:) = 0.04*ZDDU(:) +ZTSR(:) = 0.04*ZDDT(:) +ZQSR(:) = 0.04*ZDDQ(:) +ZDELTAU10N(:) = ZDDU(:) +ZDELTAT10N(:) = ZDDT(:) +ZDELTAQ10N(:) = ZDDQ(:) +JITER(:) = 99 +! +! In the following, we suppose that Richardson number PRI < XRIMAX +! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied ! +!------------------------------------------------------------------------------- +! +! 3. ITERATIVE LOOP TO COMPUTE U*, T*, Q*. +! ------------------------------------------ +! +DO JJ=1,NITERFL + DO JLON=1,SIZE(PTA) +! + IF (JCV(JLON) == -1) THEN + ZUSR0(JLON)=ZUSR(JLON) + ZTSR0(JLON)=ZTSR(JLON) + ZQSR0(JLON)=ZQSR(JLON) + IF (JJ == NITERMAX+1 .OR. JJ == NITERMAX+NITERSUP) THEN + ZDELTAU10N(JLON) = 0.5*(ZDUSTO(JLON)+ZDELTAU10N(JLON)) !forced convergence + ZDELTAT10N(JLON) = 0.5*(ZDTSTO(JLON)+ZDELTAT10N(JLON)) + ZDELTAQ10N(JLON) = 0.5*(ZDQSTO(JLON)+ZDELTAQ10N(JLON)) + IF (JJ == NITERMAX+NITERSUP) JCV(JLON)=3 + ENDIF + ZDUSTO(JLON) = ZDELTAU10N(JLON) + ZDTSTO(JLON) = ZDELTAT10N(JLON) + ZDQSTO(JLON) = ZDELTAQ10N(JLON) +! +! 3.1. Neutral parameter for wind speed (ECUME_V6 formulation) +! + IF (ZDELTAU10N(JLON) <= ZUTU) THEN + ZPARUN(JLON) = ZCOEFU(0) + ZCOEFU(1)*ZDELTAU10N(JLON) & + + ZCOEFU(2)*ZDELTAU10N(JLON)**2 & + + ZCOEFU(3)*ZDELTAU10N(JLON)**3 & + + ZCOEFU(4)*ZDELTAU10N(JLON)**4 & + + ZCOEFU(5)*ZDELTAU10N(JLON)**5 + ELSE + ZPARUN(JLON) = ZCDIRU*(ZDELTAU10N(JLON)-ZUTU) + ZORDOU + ENDIF + PCDN(JLON) = (ZPARUN(JLON)/ZDELTAU10N(JLON))**2 +! +! 3.2. Neutral parameter for temperature (ECUME_V6 formulation) +! + IF (ZDELTAU10N(JLON) <= ZUTT) THEN + ZPARTN(JLON) = ZCOEFT(0) + ZCOEFT(1)*ZDELTAU10N(JLON) & + + ZCOEFT(2)*ZDELTAU10N(JLON)**2 & + + ZCOEFT(3)*ZDELTAU10N(JLON)**3 & + + ZCOEFT(4)*ZDELTAU10N(JLON)**4 + ELSE + ZPARTN(JLON) = ZCDIRT*(ZDELTAU10N(JLON)-ZUTT) + ZORDOT + ENDIF +! +! 3.3. Neutral parameter for humidity (ECUME_V6 formulation) +! + IF (ZDELTAU10N(JLON) <= ZUTQ) THEN + ZPARQN(JLON) = ZCOEFQ(0) + ZCOEFQ(1)*ZDELTAU10N(JLON) & + + ZCOEFQ(2)*ZDELTAU10N(JLON)**2 + ELSE + ZPARQN(JLON) = ZCDIRQ*(ZDELTAU10N(JLON)-ZUTQ) + ZORDOQ + ENDIF +! +! 3.4. Scaling parameters U*, T*, Q* +! + ZUSR(JLON) = ZPARUN(JLON) + ZTSR(JLON) = ZPARTN(JLON)*ZDELTAT10N(JLON)/ZDELTAU10N(JLON) + ZQSR(JLON) = ZPARQN(JLON)*ZDELTAQ10N(JLON)/ZDELTAU10N(JLON) +! +! 3.4b Gustiness factor (Deardorff 1970) +! +! 3.4c Cool skin correction +! +! 3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979) +! +! For U + ZLMOU = PUREF(JLON)*XG*XKARMAN*(ZTSR(JLON)/PTA(JLON) & + +ZETV*ZQSR(JLON)/(1.0+ZETV*PQA(JLON)))/ZUSR(JLON)**2 +! For T/Q + ZLMOT = ZLMOU*(PZREF(JLON)/PUREF(JLON)) + ZLMOU = MAX(MIN(ZLMOU,ZLMOMAX),ZLMOMIN) + ZLMOT = MAX(MIN(ZLMOT,ZLMOMAX),ZLMOMIN) +! +! 3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970) +! Modified to include convective form following Fairall (unpublished) +! +! For U + IF (ZLMOU == 0.0) THEN + ZPSI_U = 0.0 + ELSEIF (ZLMOU > 0.0) THEN + ZPSI_U = -ZGMA*ZLMOU + ELSE + ZCHIK = (1.0-ZBTA*ZLMOU)**0.25 + ZPSIK = 2.0*LOG((1.0+ZCHIK)/2.0) & + +LOG((1.0+ZCHIK**2)/2.0) & + -2.0*ATAN(ZCHIK)+0.5*XPI + ZCHIC = (1.0-12.87*ZLMOU)**(1.0/3.0) !for very unstable conditions + ZPSIC = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0) & + -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) & + +XPI/ZSQR3 + ZPSI_U = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOU**2) !match Kansas & free-conv. forms + ENDIF + ZPSIU(JLON) = ZPSI_U +! For T/Q + IF (ZLMOT == 0.0) THEN + ZPSI_T = 0.0 + ELSEIF (ZLMOT > 0.0) THEN + ZPSI_T = -ZGMA*ZLMOT + ELSE + ZCHIK = (1.0-ZBTA*ZLMOT)**0.25 + ZPSIK = 2.0*LOG((1.0+ZCHIK**2)/2.0) + ZCHIC = (1.0-12.87*ZLMOT)**(1.0/3.0) !for very unstable conditions + ZPSIC = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0) & + -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) & + +XPI/ZSQR3 + ZPSI_T = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOT**2) !match Kansas & free-conv. forms + ENDIF + ZPSIT(JLON) = ZPSI_T +! +! 3.7. Update air-sea gradients +! + ZDDU(JLON) = ZDU(JLON) + ZDDT(JLON) = ZDT(JLON) + ZDDQ(JLON) = ZDQ(JLON) + ZDDU(JLON) = SIGN(MAX(ABS(ZDDU(JLON)),10.0*ZDUSR0),ZDDU(JLON)) + ZDDT(JLON) = SIGN(MAX(ABS(ZDDT(JLON)),10.0*ZDTSR0),ZDDT(JLON)) + ZDDQ(JLON) = SIGN(MAX(ABS(ZDDQ(JLON)),10.0*ZDQSR0),ZDDQ(JLON)) + ZLOGUS10 = LOG(PUREF(JLON)/10.0) + ZLOGTS10 = LOG(PZREF(JLON)/10.0) + ZDELTAU10N(JLON) = ZDDU(JLON)-ZUSR(JLON)*(ZLOGUS10-ZPSI_U)/XKARMAN + ZDELTAT10N(JLON) = ZDDT(JLON)-ZTSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN + ZDELTAQ10N(JLON) = ZDDQ(JLON)-ZQSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN + ZDELTAU10N(JLON) = SIGN(MAX(ABS(ZDELTAU10N(JLON)),10.0*ZDUSR0), & + ZDELTAU10N(JLON)) + ZDELTAT10N(JLON) = SIGN(MAX(ABS(ZDELTAT10N(JLON)),10.0*ZDTSR0), & + ZDELTAT10N(JLON)) + ZDELTAQ10N(JLON) = SIGN(MAX(ABS(ZDELTAQ10N(JLON)),10.0*ZDQSR0), & + ZDELTAQ10N(JLON)) +! +! 3.8. Test convergence for U*, T*, Q* +! + IF (ABS(ZUSR(JLON)-ZUSR0(JLON)) < ZDUSR0 .AND. & + ABS(ZTSR(JLON)-ZTSR0(JLON)) < ZDTSR0 .AND. & + ABS(ZQSR(JLON)-ZQSR0(JLON)) < ZDQSR0) THEN + JCV(JLON) = 1 !free convergence + IF (JJ >= NITERMAX+1) JCV(JLON) = 2 !leaded convergence + ENDIF + JITER(JLON) = JJ + ENDIF +! + ENDDO +ENDDO +! +!------------------------------------------------------------------------------- +! +! 4. COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS. +! --------------------------------------------------------------- +! +DO JLON=1,SIZE(PTA) +! +! 4.1. Surface turbulent fluxes +! (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat) +! + ZTAU(JLON) = -PRHOA(JLON)*ZUSR(JLON)**2 + ZHF(JLON) = -PRHOA(JLON)*ZCPA(JLON)*ZUSR(JLON)*ZTSR(JLON) + ZEF(JLON) = -PRHOA(JLON)*ZLVS(JLON)*ZUSR(JLON)*ZQSR(JLON) + + write(*,*) "ZTAU = ",ZTAU(JLON) + write(*,*) "SENS = ",ZHF(JLON) + write(*,*) "LAT = ",ZEF(JLON) + +! +! 4.2. Exchange coefficients PCD, PCH, PCE +! + PCD(JLON) = (ZUSR(JLON)/ZDDU(JLON))**2 + PCH(JLON) = (ZUSR(JLON)*ZTSR(JLON))/(ZDDU(JLON)*ZDDT(JLON)) + PCE(JLON) = (ZUSR(JLON)*ZQSR(JLON))/(ZDDU(JLON)*ZDDQ(JLON)) + + write(*,*) "ZUSR = ",ZUSR(JLON) + write(*,*) "ZTSR = ",ZTSR(JLON) + write(*,*) "ZQSR = ",ZQSR(JLON) + +! +! 4.3. Stochastic perturbation of turbulent fluxes +! +! IF( OPERTFLUX )THEN +! ZTAU(JLON) = ZTAU(JLON)* ( 1. + PPERTFLUX(JLON) / 2. ) +! ZHF (JLON) = ZHF(JLON)* ( 1. + PPERTFLUX(JLON) / 2. ) +! ZEF (JLON) = ZEF(JLON)* ( 1. + PPERTFLUX(JLON) / 2. ) +! ENDIF +! +ENDDO +! +!------------------------------------------------------------------------------- +! +! 5. COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL. +! (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0) +! ----------------------------------------------------- +! +IF (OPRECIP) THEN + DO JLON=1,SIZE(PTA) +! +! 5.1. Momentum flux due to rainfall (ZTAUR, N/m2) +! +! See pp3752 in FBR96. + ZTAUR(JLON) = -0.85*PRAIN(JLON)*PVMOD(JLON) +! +! 5.2. Sensible heat flux due to rainfall (ZRF, W/m2) +! +! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K), +! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity +! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB +! as wet-bulb factor (Eq.11 in GoF95). +! + ZTAC = PTA(JLON)-XTT + ZCPWA = 4217.51 -3.65566*ZTAC +0.1381*ZTAC**2 & + -2.8309E-03*ZTAC**3 +3.42061E-05*ZTAC**4 & + -2.18107E-07*ZTAC**5 +5.74535E-10*ZTAC**6 + ZDQSDT = (ZLVA(JLON)*ZQSATA(JLON))/(XRV*PTA(JLON)**2) + ZDWAT = 2.11E-05*(ZP00/PPA(JLON))*(PTA(JLON)/XTT)**1.94 + ZDTMP = (1.0+3.309E-03*ZTAC-1.44E-06*ZTAC**2) & + *0.02411/(PRHOA(JLON)*ZCPA(JLON)) + ZBULB = 1.0/(1.0+ZDQSDT*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP)) + ZRF(JLON) = PRAIN(JLON)*ZCPWA*ZBULB*((PSST(JLON)-PTA(JLON)) & + +(PQSATA(JLON)-PQA(JLON))*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP)) +! + ENDDO +ENDIF +! +!------------------------------------------------------------------------------- +! +! 6. COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2). +! ------------------------------------------------------------------------ +! +! See Eq.21 and Eq.22 in FBR96. +IF (OPWEBB) THEN + DO JLON=1,SIZE(PTA) + ZWW = (1.0+ZETV)*(ZUSR(JLON)*ZQSR(JLON)) & + +(1.0+(1.0+ZETV)*PQA(JLON))*(ZUSR(JLON)*ZTSR(JLON))/PTA(JLON) + ZEFWEBB(JLON) = -PRHOA(JLON)*ZLVS(JLON)*ZWW*PQA(JLON) + ENDDO +ENDIF +! +!------------------------------------------------------------------------------- +! +! 7. FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS. +! --------------------------------------------------------------- +! +! 7.1. Richardson number +! +! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,PQA, & +! PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI) +! +! 7.2. Friction velocity which contains correction due to rain +! +ZUSTAR2(:) = -(ZTAU(:)+ZTAUR(:))/PRHOA(:) !>>0 as ZTAU<<0 & ZTAUR<=0 +! +IF (OPRECIP) THEN + PCD(:) = ZUSTAR2(:)/ZDDU(:)**2 +ENDIF +! +PUSTAR(:) = SQRT(ZUSTAR2(:)) !>>0 +! +! 7.3. Aerodynamical conductance and resistance +! +ZAC (:) = PCH(:)*ZDDU(:) +PRESA(:) = 1.0/ZAC(:) +! +! 7.4. Total surface fluxes +! +PSFTH(:) = ZHF(:)+ZRF(:) +PSFTQ(:) = (ZEF(:)+ZEFWEBB(:))/ZLVS(:) +! +! 7.5. Charnock number +! +IF (CCHARNOCK == 'OLD') THEN + ZCHARN(:) = XVCHRNK +ELSE !modified for moderate wind speed as in COARE3.0 + ZCHARN(:) = MIN(0.018,MAX(0.011,0.011+(0.007/8.0)*(ZDDU(:)-10.0))) +ENDIF +! +! 7.6. Roughness lengths Z0 and Z0H over sea +! +!IF (KZ0 == 0) THEN ! ARPEGE formulation +! PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + XVZ0CM*PCD(:)/PCDN(:) +! PZ0HSEA(:) = PZ0SEA (:) +!ELSEIF (KZ0 == 1) THEN ! Smith (1988) formulation +! PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + 0.11*ZVISA(:)/PUSTAR(:) +! PZ0HSEA(:) = PZ0SEA (:) +!ELSEIF (KZ0 == 2) THEN ! Direct computation using the stability functions +! DO JLON=1,SIZE(PTA) +! PZ0SEA (JLON) = PUREF(JLON)/EXP(XKARMAN*ZDDU(JLON)/PUSTAR(JLON)+ZPSIU(JLON)) +! Z0TSEA = PZREF(JLON)/EXP(XKARMAN*ZDDT(JLON)/ZTSR (JLON)+ZPSIT(JLON)) +! Z0QSEA = PZREF(JLON)/EXP(XKARMAN*ZDDQ(JLON)/ZQSR (JLON)+ZPSIT(JLON)) +! PZ0HSEA(JLON) = 0.5*(Z0TSEA+Z0QSEA) +! ENDDO +!ENDIF + +write(*,*) "JLON ",JLON +write(*,*) "PTA ",klon,PTA +write(*,*) "PCD ",SIZE(PCD),PCD +write(*,*) "PCQ ",SIZE(PCE),PCE +write(*,*) "PCH ",SIZE(PCH),PCH + +coeffs = [PCD,& + PCE,& + PCH] +! +! +!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- + END SUBROUTINE ECUMEV6_FLUX Index: LMDZ6/trunk/libf/phylmd/ini_csts.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/ini_csts.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/ini_csts.F90 (revision 4722) @@ -0,0 +1,187 @@ +!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SFX_LIC for details. version 1. +! ######### + SUBROUTINE INI_CSTS +! ################## +! +!!**** *INI_CSTS * - routine to initialize the module MODD_CST +!! +!! PURPOSE +!! ------- +! The purpose of this routine is to initialize the physical constants +! stored in module MODD_CST. +! +! +!!** METHOD +!! ------ +!! The physical constants are set to their numerical values +!! +!! +!! EXTERNAL +!! -------- +!! FMLOOK : to retrieve logical unit number associated to a file +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST : contains physical constants +!! +!! REFERENCE +!! --------- +!! Book2 of the documentation (module MODD_CST, routine INI_CSTS) +!! +!! +!! AUTHOR +!! ------ +!! V. Ducrocq * Meteo France * +!! +!! MODIFICATIONS +!! ------------- +!! Original 18/05/94 +!! J. Stein 02/01/95 add the volumic mass of liquid water +!! J.-P. Pinty 13/12/95 add the water vapor pressure over solid ice +!! J. Stein 29/06/97 add XTH00 +!! V. Masson 05/10/98 add XRHOLI +!! C. Mari 31/10/00 add NDAYSEC +!! V. Masson 01/03/03 add XCONDI +!! A. Voldoire 01/12/09 add XTTSI, XICEC, XTTS for ESM +!! J. Escobar 28/03/2014 for pb with emissivity/aerosol reset XSURF_TINY=1.0e-80 in real8 case +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_CSTS +! +! +!USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +!USE PARKIND1 ,ONLY : JPRB +! +!USE MODI_INI_CTURBS +! +!USE MODI_INI_OCEAN_CSTS +! +!USE MODI_INI_SURF_CSTS +! +IMPLICIT NONE +! +!------------------------------------------------------------------------------- +! +!* 1. FUNDAMENTAL CONSTANTS +! --------------------- +! + +!REAL(KIND=JPRB) :: ZHOOK_HANDLE + +!IF (LHOOK) CALL DR_HOOK('INI_CSTS',0,ZHOOK_HANDLE) + +!#ifdef SFX_MNH +!#ifdef MNH_MPI_DOUBLE_PRECISION +XSURF_TINY = 1.0e-80 +!#else +!XSURF_TINY = TINY (XSURF_TINY ) +!#endif +!#else +XSURF_TINY = 1.0e-80 +!#endif +XSURF_TINY_12 = SQRT (XSURF_TINY ) +XSURF_EPSILON = EPSILON (XSURF_EPSILON ) * 10.0 + +XPI = 2.*ASIN(1.) +XKARMAN = 0.4 +XBOLTZ = 1.380658E-23 +XLIGHTSPEED = 299792458. +XPLANCK = 6.6260755E-34 +XAVOGADRO = 6.0221367E+23 +! +!------------------------------------------------------------------------------- +! +!* 2. ASTRONOMICAL CONSTANTS +! ---------------------- +! +XDAY = 86400. +XSIYEA = 365.25*XDAY*2.*XPI/ 6.283076 +XSIDAY = XDAY/(1.+XDAY/XSIYEA) +XOMEGA = 2.*XPI/XSIDAY +NDAYSEC = 24*3600 ! Number of seconds in a day +! +!-------------------------------------------------------------------------------! +! +! +!* 3. TERRESTRIAL GEOIDE CONSTANTS +! ---------------------------- +! +XRADIUS = 6371229. +XG = 9.80665 +! +!------------------------------------------------------------------------------- +! +!* 4. REFERENCE PRESSURE +! ------------------- +! +XP00 = 1.E5 +XTH00 = 300. +!------------------------------------------------------------------------------- +! +!* 5. RADIATION CONSTANTS +! ------------------- +! +!JUAN OVERFLOW XSTEFAN = 2.* XPI**5 * XBOLTZ**4 / (15.* XLIGHTSPEED**2 * XPLANCK**3) +XSTEFAN = ( 2.* XPI**5 / 15. ) * ( (XBOLTZ / XPLANCK)* XBOLTZ ) * (XBOLTZ/(XLIGHTSPEED*XPLANCK))**2 +XI0 = 1370. +! +!------------------------------------------------------------------------------- +! +!* 6. THERMODYNAMIC CONSTANTS +! ----------------------- +! +XMD = 28.9644E-3 +XMV = 18.0153E-3 +XRD = XAVOGADRO * XBOLTZ / XMD +XRV = XAVOGADRO * XBOLTZ / XMV +XCPD = 7.* XRD /2. +XCPV = 4.* XRV +XRHOLW = 1000. +XRHOLI = 917. +XCONDI = 2.22 +XCL = 4.218E+3 +XCI = 2.106E+3 +XTT = 273.16 +XTTSI = XTT - 1.8 +XICEC = 0.5 +XTTS = XTT*(1-XICEC) + XTTSI*XICEC +XLVTT = 2.5008E+6 +XLSTT = 2.8345E+6 +XLMTT = XLSTT - XLVTT +XESTT = 611.14 +XGAMW = (XCL - XCPV) / XRV +XBETAW = (XLVTT/XRV) + (XGAMW * XTT) +XALPW = LOG(XESTT) + (XBETAW /XTT) + (XGAMW *LOG(XTT)) +XGAMI = (XCI - XCPV) / XRV +XBETAI = (XLSTT/XRV) + (XGAMI * XTT) +XALPI = LOG(XESTT) + (XBETAI /XTT) + (XGAMI *LOG(XTT)) +! +!------------------------------------------------------------------------------- +! +!* 7. TURBULENCE CONSTANTS +! -------------------- +! +! CALL INI_CTURBS +!------------------------------------------------------------------------------- +! +!* 8. OCEAN CONSTANTS +! --------------- +! +! CALL INI_OCEAN_CSTS +! +!* 9. SURFACE CONSTANTS +! ----------------- +! +! CALL INI_SURF_CSTS +!IF (LHOOK) CALL DR_HOOK('INI_CSTS',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE INI_CSTS Index: LMDZ6/trunk/libf/phylmd/modd_csts.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/modd_csts.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/modd_csts.F90 (revision 4722) @@ -0,0 +1,91 @@ +!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SFX_LIC for details. version 1. +! ############### + MODULE MODD_CSTS +! ############### +! +!!**** *MODD_CSTS* - declaration of Physic constants +!! +!! PURPOSE +!! ------- +! The purpose of this declarative module is to declare the +! Physics constants. +! +!! +!!** IMPLICIT ARGUMENTS +!! ------------------ +!! None +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! V. Ducrocq *Meteo France* +!! +!! MODIFICATIONS +!! ------------- +!! Original 16/05/94 +!! J. Stein 02/01/95 add xrholw +!! J.-P. Pinty 13/12/95 add XALPI,XBETAI,XGAMI +!! J. Stein 25/07/97 add XTH00 +!! V. Masson 05/10/98 add XRHOLI +!! C. Mari 31/10/00 add NDAYSEC +!! J. Escobar 06/13 add XSURF_TIMY XSURF_TIMY_12 XSURF_EPSILON for REAL*4 +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +IMPLICIT NONE +REAL,SAVE :: XPI ! Pi +! +REAL,SAVE :: XDAY,XSIYEA,XSIDAY ! day duration, sideral year duration, + ! sideral day duration +! +REAL,SAVE :: XKARMAN ! von karman constant +REAL,SAVE :: XLIGHTSPEED ! light speed +REAL,SAVE :: XPLANCK ! Planck constant +REAL,SAVE :: XBOLTZ ! Boltzman constant +REAL,SAVE :: XAVOGADRO ! Avogadro number +! +REAL,SAVE :: XRADIUS,XOMEGA ! Earth radius, earth rotation +REAL,SAVE :: XG ! Gravity constant +! +REAL,SAVE :: XP00 ! Reference pressure +! +REAL,SAVE :: XSTEFAN,XI0 ! Stefan-Boltzman constant, solar constant +! +REAL,SAVE :: XMD,XMV ! Molar mass of dry air and molar mass of vapor +REAL,SAVE :: XRD,XRV ! Gaz constant for dry air, gaz constant for vapor +REAL,SAVE :: XCPD,XCPV ! Cpd (dry air), Cpv (vapor) +REAL,SAVE :: XRHOLW ! Volumic mass of liquid water +REAL,SAVE :: XCL,XCI ! Cl (liquid), Ci (ice) +REAL,SAVE :: XTT ! Triple point temperature +REAL,SAVE :: XTTSI ! Temperature of ice fusion over salty sea +REAL,SAVE :: XTTS ! Equivalent temperature of ice fusion over a mixed of sea and sea-ice +REAL,SAVE :: XICEC ! Threshold fraction over which the tile is considered as only covered with ice +REAL,SAVE :: XLVTT ! Vaporization heat constant +REAL,SAVE :: XLSTT ! Sublimation heat constant +REAL,SAVE :: XLMTT ! Melting heat constant +REAL,SAVE :: XESTT ! Saturation vapor pressure at triple point + ! temperature +REAL,SAVE :: XALPW,XBETAW,XGAMW ! Constants for saturation vapor + ! pressure function +REAL,SAVE :: XALPI,XBETAI,XGAMI ! Constants for saturation vapor + ! pressure function over solid ice +REAL, SAVE :: XTH00 ! reference value for the potential + ! temperature +REAL,SAVE :: XRHOLI ! Volumic mass of ice +REAL,SAVE :: XCONDI ! thermal conductivity of ice (W m-1 K-1) +! +INTEGER, SAVE :: NDAYSEC ! Number of seconds in a day +! +REAL,SAVE :: XSURF_TINY ! minimum real on this machine +REAL,SAVE :: XSURF_TINY_12 ! sqrt(minimum real on this machine) +REAL,SAVE :: XSURF_EPSILON ! minimum space with 1.0 +! +END MODULE MODD_CSTS + Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4722) @@ -348,5 +348,5 @@ REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle - REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall + REAL, DIMENSION(klon), INTENT(INOUT) :: rain_f ! rain fall REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall REAL, DIMENSION(klon), INTENT(IN) :: bs_f ! blowing snow fall @@ -1071,5 +1071,5 @@ !albedo SB <<< yrain_f = 0.0 ; ysnow_f = 0.0 ; ybs_f=0.0 ; yfder = 0.0 ; ysolsw = 0.0 - ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yu1 = 0.0 + ysollw = 0.0 ; yz0m = 0.0 ; yz0h = 0.0 ; yz0h_oupas = 0.0 ; yu1 = 0.0 yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0 ; yqbs1 = 0.0 ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0 @@ -1598,7 +1598,7 @@ ENDDO CALL cdrag(knon, nsrf, & - speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),& + speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, & yts, yqsurf, yz0m, yz0h, yri0, 0, & - ycdragm, ycdragh, zri1, pref ) + ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1)) ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013 @@ -1632,7 +1632,7 @@ CALL cdrag(knon, nsrf, & - speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),& + speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,& yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, & - ycdragm_x, ycdragh_x, zri1_x, pref_x ) + ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) ) ! --- special Dice. JYG+MPL 25112013 @@ -1659,7 +1659,7 @@ ENDDO CALL cdrag(knon, nsrf, & - speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),& + speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,& yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, & - ycdragm_w, ycdragh_w, zri1_w, pref_w ) + ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) ) ! !!!bug !! zgeo1(:) = wake_s(:)*zgeo1_w(:) + (1.-wake_s(:))*zgeo1_x(:) @@ -2086,5 +2086,5 @@ yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, & yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar) + yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol) ENDIF @@ -3072,5 +3072,5 @@ IF (iflag_split .eq.0) THEN IF (iflag_new_t2mq2m==1) THEN - CALL stdlevvarn(klon, knon, nsrf, zxli, & + CALL stdlevvarn(klon, knon, nsrf, zxli, & uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & @@ -3081,5 +3081,5 @@ uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m, yq2m, yt10m, yq10m, yu10m, yustar) + yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol) ENDIF ELSE !(iflag_split .eq.0) @@ -3099,9 +3099,9 @@ uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, & tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x) + yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol) CALL stdlevvar(klon, knon, nsrf, zxli, & uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, & tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, & - yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w) + yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol) ENDIF !!! Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4722) @@ -1259,4 +1259,7 @@ REAL, dimension(klon,klev) :: t_env,q_env + REAL, dimension(klon) :: pr_et + REAL, dimension(klon) :: w_et, jlr_g_c, jlr_g_s + REAL pi INTEGER ieru @@ -2804,4 +2807,15 @@ ELSE IF (iflag_gusts==2) THEN gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon) + !!!! modif olivier torres + ELSE IF (iflag_gusts==3) THEN + w_et=wstar(1,3) + jlr_g_s=(0.65*w_et)**2 + pr_et=rain_con*8640 + jlr_g_c = (((19.8*(pr_et(1:klon)**2))/(1.5+pr_et(1:klon)+pr_et(1:klon)**2))**(0.4))**2 + gustiness(1:klon)=jlr_g_c+jlr_g_s +!! write(*,*) "rain ",pr_et +!! write(*,*) "jlr_g_c",jlr_g_c +!! write(*,*) "wstar",wstar(1,3) +!! write(*,*) "jlr_g_s",jlr_g_s ! ELSE IF (iflag_gusts==2) THEN ! do i = 1, klon Index: LMDZ6/trunk/libf/phylmd/qsat_seawater.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/qsat_seawater.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/qsat_seawater.F90 (revision 4722) @@ -0,0 +1,119 @@ +!------------------------------------------------------------------------------- +! +! ###################################### +REAL FUNCTION QSAT_SEAWATER(PT,PP) +! ###################################### +! +!!**** *QSATW * - function to compute saturation vapor humidity from +!! temperature +!! +!! PURPOSE +!! ------- +! The purpose of this function is to compute the saturation vapor +! pressure from temperature over saline seawater +! +! +!!** METHOD +!! ------ +!! Given temperature T (PT), the saturation vapor pressure es(T) +!! (FOES(PT)) is computed by integration of the Clapeyron equation +!! from the triple point temperature Tt (XTT) and the saturation vapor +!! pressure of the triple point es(Tt) (XESTT), i.e +!! The reduction due to salinity is compute with the factor 0.98 (reduction of 2%) +!! +!! es(T)= 0.98*EXP( alphaw - betaw /T - gammaw Log(T) ) +!! +!! with : +!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt) +!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt +!! gammaw (XGAMW) = (Cl -Cpv) /Rv +!! +!! Then, the specific humidity at saturation is deduced. +!! +!! +!! EXTERNAL +!! -------- +!! NONE +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST : comtains physical constants +!! XALPW : Constant for saturation vapor pressure function +!! XBETAW : Constant for saturation vapor pressure function +!! XGAMW : Constant for saturation vapor pressure function +!! +!! REFERENCE +!! --------- +!! Book2 of documentation of Meso-NH +!! Zeng, X., Zhao, M., and Dickinson, R. E., 1998 : Intercomparaison of bulk +!! aerodynamic algorithm for the computation of sea surface fluxes using +!! TOGA COARE and TAO data. Journal of Climate, vol 11, n°10, pp 2628--2644 +!! +!! +!! AUTHOR +!! ------ +!! C. Lebeaupin * Meteo France * +!! +!! MODIFICATIONS +!! ------------- +!! Original 6/04/2005 +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_CSTS +USE dimphy +USE indice_sol_mod +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments and results +! +! + +REAL, DIMENSION(klon), INTENT(IN) :: PT ! Temperature + ! (Kelvin) +REAL, DIMENSION(klon), INTENT(IN) :: PP ! Pressure + ! (Pa) +REAL, DIMENSION(SIZE(PT)) :: PQSAT ! saturation vapor + ! specific humidity + ! with respect to + ! water (kg/kg) + +! +!* 0.2 Declarations of local variables +! +REAL, DIMENSION(SIZE(PT)) :: ZFOES ! saturation vapor + ! pressure + ! (Pascal) +! +REAL, DIMENSION(SIZE(PT)) :: ZWORK1 +REAL :: ZWORK2 +!REAL(KIND=JPRB) :: ZHOOK_HANDLE +!------------------------------------------------------------------------------- +! +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',0,ZHOOK_HANDLE) +! +!ZFOES = 1 !PSAT(PT(:)) +!ZFOES = 0.98*ZFOES +!ZFOES(:) = ZFOES(:)*1013.25E+02 !convert from atm to Pa +ZFOES = 0.98*EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT) ) +! vapor pressure reduction of 2% over saline seawater could have a significant +! impact on the computation of surface latent heat flux under strong wind +! conditions (Zeng et al, 1998). +! +ZWORK1 = ZFOES/PP +ZWORK2 = XRD/XRV +! +!* 2. COMPUTE SATURATION HUMIDITY +! --------------------------- +! +PQSAT = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1) +! +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',1,ZHOOK_HANDLE) +!------------------------------------------------------------------------------- +! +END FUNCTION QSAT_SEAWATER +! +!------------------------------------------------------------------------------- Index: LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90 (revision 4722) @@ -0,0 +1,104 @@ +!------------------------------------------------------------------------------- +! +! ###################################### +REAL FUNCTION QSAT_SEAWATER2(PT,PP,PSSS) +! ###################################### +! +!!**** *QSATW * - function to compute saturation vapor humidity from +!! temperature +!! +!! PURPOSE +!! ------- +! The purpose of this function is to compute the saturation vapor +! pressure from temperature over saline seawater +! +! +!!** METHOD +!! ------ +!! Given temperature T (PT) and salinity S (PSSS), the saturation vapor +!! pressure es(T,S) (FOES(PT,PSSS)) is computed following Weiss and Price +!! (1980). +!! +!! Then, the specific humidity at saturation is deduced. +!! +!! +!! EXTERNAL +!! -------- +!! NONE +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST : contains physical constants +!! +!! REFERENCE +!! --------- +!! Weiss, R.F., and Price, B.A., 1980 : Nitrous oxide solubility in water +!! and seawater. Marine Chemistry, n°8, pp 347-359. +!! +!! +!! AUTHOR +!! ------ +!! S. Belamari * Meteo France * +!! +!! MODIFICATIONS +!! ------------- +!! Original 19/03/2014 +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_CSTS, ONLY : XRD, XRV +USE dimphy +USE indice_sol_mod +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments and results +! +! +REAL, DIMENSION(klon), INTENT(IN) :: PT ! Temperature + ! (Kelvin) +REAL, DIMENSION(klon), INTENT(IN) :: PP ! Pressure + ! (Pascal) +REAL, DIMENSION(klon), INTENT(IN) :: PSSS ! Salinity + ! (g/kg) +REAL, DIMENSION(SIZE(PT)) :: PQSATA ! saturation vapor + ! specific humidity + ! with respect to + ! water (kg/kg) +! +!* 0.2 Declarations of local variables +! +REAL, DIMENSION(SIZE(PT)) :: ZFOES ! saturation vapor + ! pressure + ! (Pascal) +! +REAL, DIMENSION(SIZE(PT)) :: ZWORK1 +REAL :: ZWORK2 +!REAL(KIND=JPRB) :: ZHOOK_HANDLE +!------------------------------------------------------------------------------- +! +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW2_1D',0,ZHOOK_HANDLE) +! +!* 1. COMPUTE SATURATION VAPOR PRESSURE +! --------------------------------- +! +ZFOES(:) = EXP( 24.4543 -67.4509*(100.0/PT(:)) -4.8489*LOG(PT(:)/100.0) & + -5.44E-04*(PSSS(:)/1.00472) ) !see Sharqawy et al (2010) Eq32 p368 +ZFOES(:) = ZFOES(:)*1013.25E+02 !convert from atm to Pa +! +ZWORK1(:) = ZFOES(:)/PP(:) +ZWORK2 = XRD/XRV +! +!* 2. COMPUTE SATURATION SPECIFIC HUMIDITY +! ------------------------------------ +! +PQSATA(:) = ZWORK2*ZWORK1(:) / (1.0+(ZWORK2-1.0)*ZWORK1(:)) +! +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW2_1D',1,ZHOOK_HANDLE) +!------------------------------------------------------------------------------- +! +END FUNCTION QSAT_SEAWATER2 +! +!------------------------------------------------------------------------------- Index: LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90 (revision 4722) +++ LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90 (revision 4722) @@ -0,0 +1,113 @@ +!------------------------------------------------------------------------------- +! +! ###################################### +REAL FUNCTION QSATSEAW_1D(PT,PP) +! ###################################### +! +!!**** *QSATW * - function to compute saturation vapor humidity from +!! temperature +!! +!! PURPOSE +!! ------- +! The purpose of this function is to compute the saturation vapor +! pressure from temperature over saline seawater +! +! +!!** METHOD +!! ------ +!! Given temperature T (PT), the saturation vapor pressure es(T) +!! (FOES(PT)) is computed by integration of the Clapeyron equation +!! from the triple point temperature Tt (XTT) and the saturation vapor +!! pressure of the triple point es(Tt) (XESTT), i.e +!! The reduction due to salinity is compute with the factor 0.98 (reduction of 2%) +!! +!! es(T)= 0.98*EXP( alphaw - betaw /T - gammaw Log(T) ) +!! +!! with : +!! alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt) +!! betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt +!! gammaw (XGAMW) = (Cl -Cpv) /Rv +!! +!! Then, the specific humidity at saturation is deduced. +!! +!! +!! EXTERNAL +!! -------- +!! NONE +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST : comtains physical constants +!! XALPW : Constant for saturation vapor pressure function +!! XBETAW : Constant for saturation vapor pressure function +!! XGAMW : Constant for saturation vapor pressure function +!! +!! REFERENCE +!! --------- +!! Book2 of documentation of Meso-NH +!! Zeng, X., Zhao, M., and Dickinson, R. E., 1998 : Intercomparaison of bulk +!! aerodynamic algorithm for the computation of sea surface fluxes using +!! TOGA COARE and TAO data. Journal of Climate, vol 11, n�10, pp 2628--2644 +!! +!! +!! AUTHOR +!! ------ +!! C. Lebeaupin * Meteo France * +!! +!! MODIFICATIONS +!! ------------- +!! Original 6/04/2005 +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_CSTS +USE dimphy +USE indice_sol_mod + +IMPLICIT NONE + +!* 0.1 Declarations of arguments and results +! +! +REAL, DIMENSION(klon), INTENT(IN) :: PT ! Temperature + ! (Kelvin) +REAL, DIMENSION(klon), INTENT(IN) :: PP ! Pressure + ! (Pa) +REAL, DIMENSION(SIZE(PT)) :: PQSAT ! saturation vapor + ! specific humidity + ! with respect to + ! water (kg/kg) +! +!* 0.2 Declarations of local variables +! +REAL, DIMENSION(SIZE(PT)) :: ZFOES ! saturation vapor + ! pressure + ! (Pascal) +! +INTEGER :: JJ ! loop index +!REAL(KIND=JPRB) :: ZHOOK_HANDLE +!------------------------------------------------------------------------------- +! +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',0,ZHOOK_HANDLE) +!DO JJ = 1, SIZE(PT) +!* 1. COMPUTE SATURATION VAPOR PRESSURE +! --------------------------------- +! +ZFOES = 0.98*EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT) ) +! vapor pressure reduction of 2% over saline seawater could have a significant +! impact on the computation of surface latent heat flux under strong wind +! conditions (Zeng et al, 1998). +! +!* 2. COMPUTE SATURATION HUMIDITY +! --------------------------- +! +PQSAT = XRD/XRV*ZFOES/PP /(1.+(XRD/XRV-1.)*ZFOES/PP) +! +!ENDDO +!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',1,ZHOOK_HANDLE) +!------------------------------------------------------------------------------- +! +END FUNCTION QSATSEAW_1D +!------------------------------------ Index: LMDZ6/trunk/libf/phylmd/screenc_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/screenc_mod.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/screenc_mod.F90 (revision 4722) @@ -18,5 +18,5 @@ ts, qsurf, z0m, z0h, psol, & ustar, testar, qstar, okri, ri1, & - pref, delu, delte, delq) + pref, delu, delte, delq, s_pblh, prain, tsol, pat1) IMPLICIT NONE !----------------------------------------------------------------------- @@ -60,6 +60,13 @@ REAL, dimension(klon), intent(in) :: speed, temp, q_zref REAL, intent(in) :: zref - REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol - REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1 + REAL, dimension(klon), intent(IN) :: ts + REAL, dimension(klon), intent(in) :: qsurf, psol + REAL, dimension(klon), intent(inout):: z0m, z0h + REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1 + + REAL, dimension(klon), intent(inout) :: s_pblh + REAL, dimension(klon), intent(in) :: prain + REAL, dimension(klon), intent(in) :: tsol + REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev ! REAL, dimension(klon), intent(out) :: pref, delu, delte, delq @@ -88,7 +95,7 @@ CALL cdrag (knon, nsrf, & speed, temp, q_zref, gref, & - psol, ts, qsurf, z0m, z0h, & + psol, s_pblh, ts, qsurf, z0m, z0h, & zri_zero,0, & - cdram, cdrah, zri1, pref) + cdram, cdrah, zri1, pref, prain, tsol, pat1) DO i = 1, knon IF(ok_prescr_ust) THEN @@ -114,5 +121,5 @@ cdrm, cdrh, okri, & ri1, iri1, & - pref, delm, delh, zri1) + pref, delm, delh, zri1, s_pblh, prain, tsol, pat1) IMPLICIT NONE !----------------------------------------------------------------------- @@ -156,6 +163,11 @@ REAL, dimension(klon), intent(in) :: speed, temp, q_zref REAL, intent(in) :: zref - REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol + REAL, dimension(klon), intent(in) :: ts, qsurf, psol + REAL, dimension(klon), intent(inout) :: z0m, z0h REAL, dimension(klon), intent(in) :: cdrm, cdrh, ri1 + REAL, dimension(klon), intent(inout) :: s_pblh + REAL, dimension(klon), intent(in) :: prain + REAL, dimension(klon), intent(in) :: tsol + REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev INTEGER, INTENT(IN) :: iri1 ! Richardson de la 1ere couche ! @@ -180,7 +192,7 @@ CALL cdrag(knon, nsrf, & speed, temp, q_zref, gref, & - psol, ts, qsurf, z0m, z0h, & + psol, s_pblh, ts, qsurf, z0m, z0h, & ri1, iri1, & - cdram, cdrah, zri1, pref) + cdram, cdrah, zri1, pref, prain, tsol, pat1) DO i = 1, knon delm(i) = sqrt(cdrm(i))/sqrt(cdram(i)) Index: LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90 (revision 4720) +++ LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90 (revision 4722) @@ -19,5 +19,5 @@ u1, v1, t1, q1, z1, & ts1, qsurf, z0m, z0h, psol, pat1, & - t_2m, q_2m, t_10m, q_10m, u_10m, ustar) + t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol) IMPLICIT NONE !------------------------------------------------------------------------- @@ -61,9 +61,13 @@ LOGICAL, intent(in) :: zxli REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1 - REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h + REAL, dimension(klon), intent(in) :: qsurf + REAL, dimension(klon), intent(inout) :: z0m, z0h REAL, dimension(klon), intent(in) :: psol, pat1 ! REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m + REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh + REAL, DIMENSION(klon), INTENT(IN) :: prain + REAL, DIMENSION(klon), INTENT(IN) :: tsol !------------------------------------------------------------------------- include "flux_arp.h" @@ -120,7 +124,7 @@ CALL cdrag(knon, nsrf, & & speed, t1, q1, z1, & - & psol, ts1, qsurf, z0m, z0h, & + & psol, s_pblh, ts1, qsurf, z0m, z0h, & & zri_zero, 0, & - & cdram, cdrah, zri1, pref) + & cdram, cdrah, zri1, pref, prain, tsol, pat1) ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013 @@ -178,5 +182,5 @@ & ts1, qsurf, z0m, z0h, psol, & & ustar, testar, qstar, okri, ri1, & - & pref, delu, delte, delq) + & pref, delu, delte, delq, s_pblh ,prain, tsol, pat1) ! DO i = 1, knon @@ -280,5 +284,5 @@ & ts1, qsurf, z0m, z0h, psol, & & ustar, testar, qstar, okri, ri1, & - & pref, delu, delte, delq) + & pref, delu, delte, delq, s_pblh ,prain, tsol, pat1) ! DO i = 1, knon @@ -357,6 +361,7 @@ INTEGER, intent(in) :: klon, knon, nsrf LOGICAL, intent(in) :: zxli - REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1 - REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h + REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, ts1, z1 + REAL, dimension(klon), intent(inout) :: z0m, z0h + REAL, dimension(klon), intent(in) :: qsurf REAL, dimension(klon), intent(in) :: psol, pat1 ! @@ -371,4 +376,7 @@ REAL, dimension(klon) :: cdmn2m, cdhn2m, fm2m, fh2m REAL, dimension(klon) :: ri2m_new + REAL, DIMENSION(klon) :: s_pblh + REAL, DIMENSION(klon) :: prain + REAL, DIMENSION(klon) :: tsol !------------------------------------------------------------------------- include "flux_arp.h" @@ -444,7 +452,7 @@ CALL cdrag(knon, nsrf, & & speed, t1, q1, z1, & - & psol, ts1, qsurf, z0m, z0h, & + & psol, s_pblh, ts1, qsurf, z0m, z0h, & & zri_zero, 0, & - & cdram, cdrah, zri1, pref) + & cdram, cdrah, zri1, pref, prain, tsol, pat1) ! @@ -468,5 +476,6 @@ & cdram, cdrah, okri, & & ri1, 1, & - & pref_new, delm_new, delh_new, ri2m) + & pref_new, delm_new, delh_new, ri2m, & + & s_pblh, prain, tsol, pat1 ) ! DO i = 1, knon @@ -535,5 +544,6 @@ & cdram, cdrah, okri, & & ri1, 0, & - & pref, delm, delh, ri2m) + & pref, delm, delh, ri2m, & + & s_pblh, prain, tsol, pat1 ) ! DO i = 1, knon @@ -614,5 +624,6 @@ & cdram, cdrah, okri, & & ri1, 1, & - & pref_new, delm_new, delh_new, ri10m) + & pref_new, delm_new, delh_new, ri10m, & + & s_pblh, prain, tsol, pat1 ) ! DO i = 1, knon @@ -671,5 +682,6 @@ & cdram, cdrah, okri, & & ri1, 0, & - & pref, delm, delh, ri10m) + & pref, delm, delh, ri10m, & + & s_pblh, prain, tsol, pat1 ) ! DO i = 1, knon