Index: /LMDZ6/trunk/arch/arch-X64_ADASTRA-GNU.fcm =================================================================== --- /LMDZ6/trunk/arch/arch-X64_ADASTRA-GNU.fcm (revision 5060) +++ /LMDZ6/trunk/arch/arch-X64_ADASTRA-GNU.fcm (revision 5061) @@ -9,5 +9,5 @@ %FPP_DEF NC_DOUBLE -%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch -march=native -fPIC +%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch -fimplicit-none -march=native -fPIC %BASE_CFLAGS -w -std=c++11 -D__XIOS_EXCEPTION # xios # /!\ LD must be written in Makefile syntax @@ -30,3 +30,2 @@ %CPP cpp # xios -%MAKE make # xios Index: /LMDZ6/trunk/arch/arch-local-gfortran-parallel.fcm =================================================================== --- /LMDZ6/trunk/arch/arch-local-gfortran-parallel.fcm (revision 5060) +++ /LMDZ6/trunk/arch/arch-local-gfortran-parallel.fcm (revision 5061) @@ -9,5 +9,5 @@ %FPP_DEF NC_DOUBLE -%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch +%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch -fimplicit-none %BASE_CFLAGS -w -std=c++11 -D__XIOS_EXCEPTION # xios # /!\ LD must be written in Makefile syntax @@ -30,5 +30,4 @@ %CPP cpp # xios -%MAKE make # xios Index: /LMDZ6/trunk/arch/arch-local-gfortran.fcm =================================================================== --- /LMDZ6/trunk/arch/arch-local-gfortran.fcm (revision 5060) +++ /LMDZ6/trunk/arch/arch-local-gfortran.fcm (revision 5061) @@ -7,5 +7,5 @@ %FPP_FLAGS -P -traditional %FPP_DEF NC_DOUBLE -%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch +%BASE_FFLAGS -cpp -ffree-line-length-0 -fdefault-real-8 -DNC_DOUBLE -fallow-argument-mismatch -fimplicit-none %PROD_FFLAGS -O3 -march=native %DEV_FFLAGS -Wall -fbounds-check Index: /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm_mod.F90 (revision 5061) +++ /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm_mod.F90 (revision 5061) @@ -0,0 +1,598 @@ +!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 coare30_flux_cnrm_mod + implicit none + private + public COARE30_FLUX_CNRM + +contains + + +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 coare_cp_mod, ONLY: PSIFCTT=>psit_30, PSIFCTU=>psiuo + +!-------------------------------- + +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 +! +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(1)/ZL10(1))) +ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1))) +ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1))) +! +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>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(1)) + ZPTZ(:) = PSIFCTT(ZZTL(1)) + ! + 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 + +end module coare30_flux_cnrm_mod Index: /LMDZ6/trunk/libf/phylmd/coare_cp_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/coare_cp_mod.F90 (revision 5061) +++ /LMDZ6/trunk/libf/phylmd/coare_cp_mod.F90 (revision 5061) @@ -0,0 +1,249 @@ +module coare_cp_mod + implicit none + private + public psit_30, psiuo, coare_cp + +contains + + real function psit_30(zet) + implicit none + real, intent(in) :: zet + + real :: x, psik, psic, f, c + + 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) + implicit none + real, intent(in) :: zet + + real :: x, psik, psic, f, c + + 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 + + +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 :: 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 < 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 > 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 > 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 + +end module coare_cp_mod