Changeset 5060


Ignore:
Timestamp:
Jul 16, 2024, 4:57:48 PM (8 weeks ago)
Author:
abarral
Message:

Refactor coare_cp.F90 & coare30_flux_cnrm.F90 into modules to reduce duplicated code. Accordingly adapt casting for coare30_flux_cnrm calls, which were inconsistent with declared interfaces.
Update deprecated math operators.
Lint end-of-line trailing spaces.
Remove dangling get_var1 from dynredem_mod.F90.

Location:
LMDZ6/trunk/libf
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/dynredem_mod.F90

    r2299 r5060  
    44  PRIVATE
    55  PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err
    6   PUBLIC :: cre_var, get_var1, put_var1, put_var2, fil, modname, msg
     6  PUBLIC :: cre_var, put_var1, put_var2, fil, modname, msg
    77  include "dimensions.h"
    88  include "paramet.h"

  • LMDZ6/trunk/libf/dyn3dmem/dynredem_mod.F90

    r2299 r5060  
    77  PRIVATE
    88  PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err
    9   PUBLIC :: cre_var, get_var1, put_var, fil, modname, msg
     9  PUBLIC :: cre_var, put_var, fil, modname, msg
    1010  CHARACTER(LEN=256), SAVE :: fil, modname
    1111  INTEGER,            SAVE :: nvarid

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

    r4777 r5060  
    2323
    2424  USE dimphy
     25  USE coare_cp_mod, ONLY: coare_cp
     26  USE coare30_flux_cnrm_mod, ONLY: coare30_flux_cnrm
    2527  USE indice_sol_mod
    2628  USE print_control_mod, ONLY: lunout, prt_level
    … …  
    341343         LPWG    = .false.
    342344         call ini_csts
    343          call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i),  &
    344              sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, &
    345              PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
    346              PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs)
     345         block
     346           real, dimension(1) :: z0m_1d, z_0h_1d, sqrt_zdu2_1d, zgeop1_rg_1d  ! convert scalar to 1D for call
     347           z0m_1d = z0m
     348           z_0h_1d = z0h
     349           sqrt_zdu2_1d = sqrt(zdu2)
     350           zgeop1_rg_1d=zgeop1(i)/RG
     351           call coare30_flux_cnrm(z0m_1d,t1(i),tsurf(i), q1(i),  &
     352               sqrt_zdu2_1d,zgeop1_rg_1d,zgeop1_rg_1d,psol(i),qsurf(i),PQSAT, &
     353               PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
     354               PRESA,prain,pat1(i),z_0h_1d, LPRECIP, LPWG, coeffs)
     355
     356         end block
    347357         cdmm(i) = coeffs(1)
    348358         cdhh(i) = coeffs(2)

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

    r4722 r5060  
    55!     #########
    66
    7 !------------------------ Modif Olivier Torres pour rajout fonction et pas appel ----------------
    8 
    9 
    10 
    11 real function PSIFCTT(zet)
    12 
    13   real, intent(in) :: zet
    14 
    15   if(zet<0) then
    16      x=(1.-(15*zet))**.5
    17      psik=2*log((1+x)/2)
    18      x=(1.-(34.15*zet))**.3333
    19      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    20      f=zet*zet/(1+zet*zet)
    21      PSIFCTT=(1-f)*psik+f*psic   
    22 
    23   else
    24      c=min(50.,.35*zet)
    25      PSIFCTT=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
    26   endif
    27 end FUNCTION PSIFCTT
    28 
    29 real function PSIFCTU(zet)
    30 
    31   real, intent(in) :: zet
    32 
    33   if (zet<0) then
    34 
    35      x=(1.-15.*zet)**.25
    36      psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
    37      x=(1.-10.15*zet)**.3333
    38      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    39      f=zet*zet/(1+zet*zet)
    40      PSIFCTU=(1-f)*psik+f*psic                                               
    41   else
    42      c=min(50.,.35*zet)
    43      PSIFCTU=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
    44   endif
    45 END FUNCTION PSIFCTU
    46 
    47 !----------------------------------- Fin Modif ---------------------------------------------------
     7module coare30_flux_cnrm_mod
     8  implicit none
     9  private
     10  public COARE30_FLUX_CNRM
     11
     12contains
    4813
    4914
    5015SUBROUTINE COARE30_FLUX_CNRM(PZ0SEA,PTA,PSST,PQA,  &
    5116            PVMOD,PZREF,PUREF,PPS,PQSATA,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,&
    52             PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs ) 
     17            PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs )
    5318!     #######################################################################
    5419!
    5520!
    56 !!****  *COARE25_FLUX* 
     21!!****  *COARE25_FLUX*
    5722!!
    5823!!    PURPOSE
    5924!!    -------
    6025!      Calculate the surface fluxes of heat, moisture, and momentum over
    61 !      sea surface with bulk algorithm COARE3.0. 
    62 !     
     26!      sea surface with bulk algorithm COARE3.0.
     27!
    6328!!**  METHOD
    6429!!    ------
    6530!      transfer coefficients were obtained using a dataset which combined COARE
    6631!      data with those from three other ETL field experiments, and reanalysis of
    67 !      the HEXMAX data (DeCosmos et al. 1996). 
    68 !      ITERMAX=3 
    69 !      Take account of the surface gravity waves on the velocity roughness and 
     32!      the HEXMAX data (DeCosmos et al. 1996).
     33!      ITERMAX=3
     34!      Take account of the surface gravity waves on the velocity roughness and
    7035!      hence the momentum transfer coefficient
    7136!        NGRVWAVES=0 no gravity waves action (Charnock) !default value
    … …  
    7742!!
    7843!!    IMPLICIT ARGUMENTS
    79 !!    ------------------ 
    80 !!     
     44!!    ------------------
     45!!
    8146!!    REFERENCE
    8247!!    ---------
    … …  
    8550!!      Gosnell et al (1995), JGR, 437-442
    8651!!      Fairall et al (1996), JGR, 1295-1308
    87 !!     
     52!!
    8853!!    AUTHOR
    8954!!    ------
    … …  
    9560!!      B. Decharme    06/2009 limitation of Ri
    9661!!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
    97 !!      B. Decharme    06/2013 bug in z0 (output) computation 
     62!!      B. Decharme    06/2013 bug in z0 (output) computation
    9863!!      J.Escobar      06/2013  for REAL4/8 add EPSILON management
    9964!!      C. Lebeaupin   03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value
    … …  
    11176USE dimphy
    11277USE indice_sol_mod
     78USE coare_cp_mod, ONLY: PSIFCTT=>psit_30, PSIFCTU=>psiuo
    11379
    11480!--------------------------------
    … …  
    159125!
    160126REAL, DIMENSION(klon), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
    161 !                                                                                 
     127!
    162128!  surface fluxes : latent heat, sensible heat, friction fluxes
    163129REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
    … …  
    176142
    177143LOGICAL,            INTENT(IN)    :: LPRECIP   !
    178 LOGICAL,            INTENT(IN)    :: LPWG    ! 
     144LOGICAL,            INTENT(IN)    :: LPWG    !
    179145real, dimension(3), intent(inout)      :: coeffs
    180146!
    … …  
    195161REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
    196162!
    197 REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref 
     163REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref
    198164REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
    199165!
    … …  
    204170REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
    205171!
    206 REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height) 
     172REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height)
    207173REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
    208174REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
    … …  
    215181REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
    216182!
    217 REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability 
     183REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability
    218184                                                     !param. z/l for u,T,q
    219185REAL, DIMENSION(SIZE(PTA))      :: ZRR
    220186REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
    221 REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q 
     187REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q
    222188!
    223189REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
    … …  
    235201REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
    236202REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
    237 REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature 
     203REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature
    238204!
    239205REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
    … …  
    254220REAL    :: QSAT_SEAWATER
    255221REAL    :: QSATSEAW_1D
    256 REAL, EXTERNAL :: PSIFCTU, PSIFCTT
    257222!
    258223INTEGER :: J, JLOOP    !loop indice
    … …  
    264229
    265230REAL       :: XVCHRNK = 0.021
    266 REAL       :: XVZ0CM = 1.0E-5 
     231REAL       :: XVZ0CM = 1.0E-5
    267232!REAL       :: XRIMAX
    268233
    … …  
    320285!
    321286!-------------------------------------------------------------------------------
    322 !       2. INITIAL GUESS FOR THE ITERATIVE METHOD 
     287!       2. INITIAL GUESS FOR THE ITERATIVE METHOD
    323288!          -------------------------------------
    324289!
    325 !       2.0     Temperature 
     290!       2.0     Temperature
    326291!
    327292! Set a non-zero value for the temperature gradient
    328293!
    329 WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.) 
     294WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.)
    330295      ZTA(:)=PTA(:)-1E-3
    331296ELSEWHERE
    332       ZTA(:)=PTA(:)     
     297      ZTA(:)=PTA(:)
    333298ENDWHERE
    334299
    335 !       2.1     Wind and humidity 
    336 !
    337 ! Sea surface specific humidity 
    338 !
    339 !PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:)) 
    340 PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:)) 
    341 
    342 
    343        
    344 !             
    345 ! Set a minimum value to wind 
     300!       2.1     Wind and humidity
     301!
     302! Sea surface specific humidity
     303!
     304!PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:))
     305PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:))
     306
     307
     308
     309!
     310! Set a minimum value to wind
    346311!
    347312!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))
    … …  
    352317
    353318!
    354 ! Specific humidity at saturation at the atm. level 
     319! Specific humidity at saturation at the atm. level
    355320!
    356321ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
    357 !ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:)) 
    358 ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:)) 
     322!ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:))
     323ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:))
    359324
    360325
    … …  
    365330IF (LPWG) ZWG(:) = 0.5
    366331!
    367 ZCHARN(:) = 0.011 
     332ZCHARN(:) = 0.011
    368333!
    369334DO J=1,SIZE(PTA)
    370335  !
    371336  !      2.2       initial guess
    372   !   
     337  !
    373338  ZDU(J) = ZVMOD(J)   !wind speed difference with surface current(=0) (m/s)
    374339                      !initial guess for gustiness factor
    … …  
    384349  ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+&
    385350             8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11
    386   ! 
     351  !
    387352  ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
    388353  ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
    … …  
    400365  !     &((ZTA(J)-XTT)*ZDUWG(J)**2)
    401366  ZRIBU(J)  = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/&
    402                (ZTA(J)*ZDUWG(J)**2) 
     367               (ZTA(J)*ZDUWG(J)**2)
    403368  !
    404369  IF (ZRIBU(J)<0.) THEN
    405370    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J))    !Unstable G and F
    406371  ELSE
    407     ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable 
     372    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable
    408373  ENDIF
    409374  !
    … …  
    413378!
    414379!  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
    415 ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF/ZL10))
    416 ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10))
    417 ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10))
     380ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(1)/ZL10(1)))
     381ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
     382ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
    418383!
    419384ZZL(:) = 0.0
    … …  
    431396  IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
    432397  !
    433   !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 
     398  !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR
    434399  !                -------------------------------------------
    435400  !
    … …  
    441406ENDDO
    442407!
    443    
     408
    444409!
    445410DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
    … …  
    447412  DO J=1,SIZE(PTA)
    448413    !
    449     IF (JLOOP.GT.ITERMAX(J)) CYCLE
     414    IF (JLOOP>ITERMAX(J)) CYCLE
    450415    !
    451416    IF (NGRVWAVES==0) THEN
    … …  
    453418    ELSE IF (NGRVWAVES==1) THEN
    454419      ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 &
    455               + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002 
     420              + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002
    456421    ELSE IF (NGRVWAVES==2) THEN
    457422      ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 &
    458               + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001 
     423              + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001
    459424    ENDIF
    460425    !
    … …  
    465430    ZZL(J) = XKARMAN * XG * PUREF(J) * &
    466431              ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / &
    467               ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) ) 
    468     ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T 
     432              ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) )
     433    ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T
    469434!    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
    470435  ENDDO
    471436  !
    472   ZPUZ(:) = PSIFCTU(ZZL(:))     
    473   ZPTZ(:) = PSIFCTT(ZZTL(:))
     437  ZPUZ(:) = PSIFCTU(ZZL(1))
     438  ZPTZ(:) = PSIFCTT(ZZTL(1))
    474439  !
    475440  DO J=1,SIZE(PTA)
    476441    !
    477     ! ZPQZ(J)=PSIFCTT(ZZQL(J))   
     442    ! ZPQZ(J)=PSIFCTT(ZZQL(J))
    478443    ZPQZ(J) = ZPTZ(J)
    479444    !
    … …  
    493458        ZWG(J) = 0.2
    494459      ENDIF
    495     ENDIF 
     460    ENDIF
    496461    ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2)
    497462    !
    … …  
    515480  !
    516481  !
    517   !            4. transfert coefficients PCD, PCH and PCE 
    518   !                 and neutral PCDN, ZCHN, ZCEN 
     482  !            4. transfert coefficients PCD, PCH and PCE
     483  !                 and neutral PCDN, ZCHN, ZCEN
    519484  !
    520485  PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
    … …  
    528493  ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT)
    529494  !
    530   !            4. 2 surface fluxes 
    531   !
    532   IF (ABS(PCDN(J))>1.E-2) THEN   !!!! secure COARE3.0 CODE 
     495  !            4. 2 surface fluxes
     496  !
     497  IF (ABS(PCDN(J))>1.E-2) THEN   !!!! secure COARE3.0 CODE
    533498    write(*,*) 'pb PCDN in COARE30: ',PCDN(J)
    534499    write(*,*) 'point: ',J,"/",SIZE(PTA)
    … …  
    543508    ZHF(J)  =  PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J)
    544509    ZEF(J)  =  PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J)
    545     !   
     510    !
    546511    !           4.3 Contributions to surface  fluxes due to rainfall
    547512    !
    … …  
    550515    !     sec) au cas humide.
    551516    IF (LPRECIP) THEN
    552       ! 
     517      !
    553518      ! heat surface  fluxes
    554519      !
    … …  
    561526      !
    562527      ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94           ! water vapour diffusivity from eq (13.3)
    563       !                                                                 ! of Pruppacher and Klett (1978)     
     528      !                                                                 ! of Pruppacher and Klett (1978)
    564529      ZALFAC(J)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
    565                    ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)     
     530                   ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)
    566531      ZCPLW(J) = 4224.8482 + ZTAC(J) * &
    567532                              ( -4.707 + ZTAC(J) * &
    568533                                (0.08499 + ZTAC(J) * &
    569534                                  (1.2826e-3 + ZTAC(J) * &
    570                                     (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat 
    571       !       
     535                                    (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat
     536      !
    572537      ZRF(J)   = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * &                    !Eq.12 in GoF95 !SIGNE?
    573538                   (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD )
    574539      !
    575       ! Momentum flux due to rainfall 
     540      ! Momentum flux due to rainfall
    576541      !
    577542      ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96
    … …  
    580545    !
    581546    !             4.4   Webb correction to latent heat flux
    582     ! 
     547    !
    583548    ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) &
    584                - ZUSR(J)*ZTSR(J)/ZTA(J)                        ! Eq.21*rhoa in FBR96   
    585     !
    586     !             4.5   friction velocity which contains correction du to rain           
     549               - ZUSR(J)*ZTSR(J)/ZTA(J)                        ! Eq.21*rhoa in FBR96
     550    !
     551    !             4.5   friction velocity which contains correction du to rain
    587552    !
    588553    ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J)
    … …  
    590555    !
    591556    !             4.6   Total surface fluxes
    592     !           
     557    !
    593558    PSFTH (J) =  ZHF(J) + ZRF(J)
    594559    PSFTQ (J) =  ZEF(J) / ZLV(J)
    595     ! 
     560    !
    596561  ENDIF
    597 ENDDO 
     562ENDDO
    598563
    599564
    … …  
    601566       PCE,&
    602567       PCH]
    603                      
     568
    604569!-------------------------------------------------------------------------------
    605570!
    606 !       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 
     571!       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS
    607572!           -----------
    608573!       5.1    Richardson number
    609 !             
     574!
    610575!
    611576!------------STOP LA --------------------
    612577!ZDIRCOSZW(:) = 1.
    613578! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,&
    614 !                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI   ) 
     579!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI   )
    615580!!
    616581!!       5.2     Aerodynamical conductance and resistance
    617 !!             
     582!!
    618583!ZAC(:) = PCH(:)*ZVMOD(:)
    619584!PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON)
    … …  
    630595!
    631596END SUBROUTINE COARE30_FLUX_CNRM
     597
     598end module coare30_flux_cnrm_mod

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

    r4722 r5060  
    1 real function psit_30(zet)
    2 
    3   real, intent(in) :: zet
    4 
    5   if(zet<0) then
    6      x=(1.-(15*zet))**.5
    7      psik=2*log((1+x)/2)
    8      x=(1.-(34.15*zet))**.3333
    9      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    10      f=zet*zet/(1+zet*zet)
    11      psit_30=(1-f)*psik+f*psic   
    12 
    13   else
    14      c=min(50.,.35*zet)
    15      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
     1module coare_cp_mod
     2  implicit none
     3  private
     4  public psit_30, psiuo, coare_cp
     5
     6contains
     7
     8  real function psit_30(zet)
     9    implicit none
     10    real, intent(in) :: zet
     11
     12    real :: x, psik, psic, f, c
     13
     14    if(zet<0) then
     15      x=(1.-(15.*zet))**.5
     16      psik=2.*log((1.+x)/2)
     17      x=(1.-(34.15*zet))**.3333
     18      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
     19      f=zet*zet/(1.+zet*zet)
     20      psit_30=(1.-f)*psik+f*psic
     21    else
     22      c=min(50.,.35*zet)
     23      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
     24    endif
     25
     26  end function psit_30
     27
     28  real function psiuo(zet)
     29    implicit none
     30    real, intent(in) :: zet
     31
     32    real :: x, psik, psic, f, c
     33
     34    if (zet<0) then
     35       x=(1.-15.*zet)**.25
     36       psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
     37       x=(1.-10.15*zet)**.3333
     38       psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
     39       f=zet*zet/(1+zet*zet)
     40       psiuo=(1-f)*psik+f*psic
     41    else
     42     c=min(50.,.35*zet)
     43     psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
    1644  endif
    17 end FUNCTION psit_30
    18 
    19 real function psiuo(zet)
    20 
    21   real, intent(in) :: zet
    22 
    23   if (zet<0) then
    24 
    25      x=(1.-15.*zet)**.25
    26      psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
    27      x=(1.-10.15*zet)**.3333
    28      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    29      f=zet*zet/(1+zet*zet)
    30      psiuo=(1-f)*psik+f*psic                                               
    31   else
    32      c=min(50.,.35*zet)
    33      psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
    34   endif
    35 END FUNCTION psiuo
    36 
    37 
    38 real function gen_q1(ts,p,dq)
    39 
    40   real, intent(in) :: ts,p,dq
    41   real es,qs
    42 
    43   es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p)
    44   qs = es*62.197/(p-37.8*es)
    45 
    46   q1 = dq + qs 
    47 
    48 end function gen_q1
    49 
    50 
     45
     46end function psiuo
    5147
    5248
    … …  
    6157  !version with shortened iteration  modified Rt and Rq
    6258  !uses wave information wave period in s and wave ht in m
    63   !no wave, standard coare 2.6 charnock:  jwave=0 
     59  !no wave, standard coare 2.6 charnock:  jwave=0
    6460  !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
    6561  !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2
    … …  
    104100  real old_usr, old_tsr, old_qsr,tmp
    105101
    106   real, external :: psit_30, psiuo, grv
     102  real, external :: grv
    107103  integer i,j,k
    108104!---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
    … …  
    113109!---------------------------------------------------------------------------------------------
    114110
    115   Ribcu=-zu/zi/.004/Beta**3 
    116   cpa=1004.67 
     111  Ribcu=-zu/zi/.004/Beta**3
     112  cpa=1004.67
    117113
    118114  t_c = t - 273.15
    119115
    120116
    121   Le=(2.501-.00237*(t_c-dt))*1e6 
    122 
    123 
    124   visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 
    125 
    126 
    127   cpv=cpa*(1+0.84*Q) 
    128   rhoa=P/(Rgas*t*(1+0.61*Q)) 
    129 
    130 
    131 
    132 
    133 
    134 
    135   ug=.2 
    136 
    137   ut=sqrt(du*du+ug*ug) 
     117  Le=(2.501-.00237*(t_c-dt))*1e6
     118
     119
     120  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3)
     121
     122
     123  cpv=cpa*(1+0.84*Q)
     124  rhoa=P/(Rgas*t*(1+0.61*Q))
     125
     126
     127
     128
     129
     130
     131  ug=.2
     132
     133  ut=sqrt(du*du+ug*ug)
    138134  ut=MAX(ut , 0.1 * MIN(10.,zu) )
    139135
    140   u10=ut*log(10/1e-4)/log(zu/1e-4) 
     136  u10=ut*log(10/1e-4)/log(zu/1e-4)
    141137  usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
    142138  zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
    143   Cd10=(von/log(10/zom10))**2 
     139  Cd10=(von/log(10/zom10))**2
    144140!  zoh10=0.40*visa/usr+1.4e-5
    145 !  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 
    146   Ch10=0.00115 
    147   Ct10=Ch10/sqrt(Cd10) 
     141!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca
     142  Ch10=0.00115
     143  Ct10=Ch10/sqrt(Cd10)
    148144  zot10=10/exp(von/Ct10)                    ! roughness length for t
    149   Cd=(von/log(zu/zom10))**2 
    150   Ct=von/log(zt/zot10) 
    151   CC=von*Ct/Cd 
     145  Cd=(von/log(zu/zom10))**2
     146  Ct=von/log(zt/zot10)
     147  CC=von*Ct/Cd
    152148
    153149  ut0 = ut
    … …  
    155151
    156152  ut = ut0
    157   Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 
    158 
    159   if (Ribu .LT. 0) then
    160      zetu=CC*Ribu/(1+Ribu/Ribcu) 
    161   else 
    162      zetu=CC*Ribu*(1+27/9*Ribu/CC) 
     153  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2
     154
     155  if (Ribu < 0) then
     156     zetu=CC*Ribu/(1+Ribu/Ribcu)
     157  else
     158     zetu=CC*Ribu*(1+27/9*Ribu/CC)
    163159  endif
    164160
    165   L10=zu/zetu 
    166 
    167   ! if (zetu .GT. 50) then 
    168   !    nits=1 
     161  L10=zu/zetu
     162
     163  ! if (zetu .GT. 50) then
     164  !    nits=1
    169165  ! endif
    170166
    171167  usr=ut*von/(log(zu/zom10)-psiuo(zetu))
    172   tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 
    173   qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 
     168  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10))
     169  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10))
    174170  !  tkt=.001
    175171
    176172  ! charnock constant - lin par morceau - constant
    177173  if ( ut <= 10. ) then
    178      charn=0.011 
     174     charn=0.011
    179175  else
    180      if (ut .GT. 18) then
     176     if (ut > 18) then
    181177        charn=0.018
    182178     else
    183         charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 
     179        charn=0.011+(ut-10)/(18-10)*(0.018-0.011)
    184180     endif
    185181  endif
    … …  
    192188
    193189  !***************  bulk loop ************
    194   do i=1, nits 
     190  do i=1, nits
    195191
    196192     zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)
    … …  
    200196!     ELSE IF (NGRVWAVES==1) THEN
    201197!       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
    202 !               + 0.11*visa/usr                       !Oost et al. 2002 
     198!               + 0.11*visa/usr                       !Oost et al. 2002
    203199!     ELSE IF (NGRVWAVES==2) THEN
    204200!       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
    205 !               + 0.11*visa/usr                       !Taulor and Yelland 2001 
    206 !     ENDIF 
    207 
    208      zom=charn*usr*usr/grav+0.11*visa/usr 
    209 
    210      rr=zom*usr/visa 
    211      L=zu/zet 
     201!               + 0.11*visa/usr                       !Taulor and Yelland 2001
     202!     ENDIF
     203
     204     zom=charn*usr*usr/grav+0.11*visa/usr
     205
     206     rr=zom*usr/visa
     207     L=zu/zet
    212208     zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
    213209     zot=zoq                         ! a modifier
    214 !     zot=0.40*visa/usr+1.4e-5 
     210!     zot=0.40*visa/usr+1.4e-5
    215211
    216212     old_usr = usr
    … …  
    218214     old_qsr = qsr
    219215
    220      usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 
    221      tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 
    222      qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 
    223 
    224      Bf=-grav/t*usr*(tsr+.61*t*qsr) 
    225 
    226      if (Bf .GT. 0) then
    227         ug=Beta*(Bf*zi)**.333 
    228      else 
    229         ug=.2 
     216     usr=ut*von/(log(zu/zom)-psiuo(zu/L))
     217     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L))
     218     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L))
     219
     220     Bf=-grav/t*usr*(tsr+.61*t*qsr)
     221
     222     if (Bf > 0) then
     223        ug=Beta*(Bf*zi)**.333
     224     else
     225        ug=.2
    230226     endif
    231227
    232      ut=sqrt(du*du+ug*ug) 
     228     ut=sqrt(du*du+ug*ug)
    233229     ut=MAX(ut , 0.1 * MIN(10.,zu) )
    234230
    … …  
    237233
    238234  ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
    239   ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 
    240   ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 
     235  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr
     236  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr
    241237
    242238  tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du
    … …  
    249245  rugosh = zot
    250246
    251 
    252247end subroutine coare_cp
     248
     249end module coare_cp_mod
Note: See TracChangeset for help on using the changeset viewer.