Ignore:
Timestamp:
Jul 22, 2024, 9:29:09 PM (4 months ago)
Author:
abarral
Message:

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecumev6_flux.F90

    r4722 r5099  
    8080!!!
    8181!-------------------------------------------------------------------------------
    82 !
     82
    8383!       0.   DECLARATIONS
    8484!            ------------
    85 !
     85
    8686USE dimphy
    8787USE indice_sol_mod
    … …  
    9393!USE MODD_SURF_ATM,         ONLY : XVCHRNK, XVZ0CM
    9494!USE MODD_REPROD_OPER,      ONLY : CCHARNOCK
    95 !
     95
    9696!USE MODE_THERMOS
    9797!USE MODI_WIND_THRESHOLD
    9898!USE MODI_SURFACE_RI
    99 !
     99
    100100!USE YOMHOOK,   ONLY : LHOOK,   DR_HOOK
    101101!USE PARKIND1,  ONLY : JPRB
    102 !
     102
    103103!USE MODI_ABOR1_SFX
    104 !
     104
    105105IMPLICIT NONE
    106 !
     106
    107107!       0.1. Declarations of arguments
    108 !
     108
    109109REAL, DIMENSION(klon), INTENT(IN)    :: PVMOD      ! module of wind at atm level (m/s)
    110110REAL, DIMENSION(klon), INTENT(IN)    :: PTA        ! air temperature at atm level (K)
    … …  
    126126LOGICAL,            INTENT(IN)    :: OPERTFLUX
    127127REAL, DIMENSION(klon), INTENT(IN)    :: PRAIN     ! precipitation rate (kg/s/m2)
    128 !
     128
    129129!INTEGER,            INTENT(IN)    :: KZ0
    130 !
     130
    131131REAL, DIMENSION(klon), INTENT(INOUT) :: PSST       ! Sea Surface Temperature (K)
    132132REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA     ! roughness length over sea
    … …  
    148148
    149149!       0.2. Declarations of local variables
    150 !
     150
    151151! specif SB
    152152INTEGER, DIMENSION(SIZE(PTA))     :: JCV        ! convergence index
    … …  
    155155REAL, DIMENSION(SIZE(PTA))        :: PEXNA      ! Exner function at atm level
    156156REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
    157 !
     157
    158158REAL, DIMENSION(SIZE(PTA))        :: ZTAU       ! momentum flux (N/m2)
    159159REAL, DIMENSION(SIZE(PTA))        :: ZHF        ! sensible heat flux (W/m2)
    … …  
    249249
    250250!REAL(KIND=JPRB) :: ZHOOK_HANDLE
    251 !
     251
    252252!-------------------------------------------------------------------------------
    253253!----------------------- Modif Olive calcul de PRHOA ---------------------------
    … …  
    280280
    281281!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',0,ZHOOK_HANDLE)
    282 !
     282
    283283ZDUSR0   = 1.E-06
    284284ZDTSR0   = 1.E-06
    285285ZDQSR0   = 1.E-09
    286 !
     286
    287287NITERMAX = 5
    288288NITERSUP = 5
    289289OPCVFLX  = .TRUE.
    290 !
     290
    291291NITERFL = NITERMAX
    292292IF (OPCVFLX) NITERFL = NITERMAX+NITERSUP
    293 !
     293
    294294ZCOEFU = (/ 1.00E-03, 3.66E-02, -1.92E-03, 2.32E-04, -7.02E-06,  6.40E-08 /)
    295295ZCOEFT = (/ 5.36E-03, 2.90E-02, -1.24E-03, 4.50E-04, -2.06E-05,       0.0 /)
    296296ZCOEFQ = (/ 1.00E-03, 3.59E-02, -2.87E-04,      0.0,       0.0,       0.0 /)
    297 !
     297
    298298ZUTU = 40.0
    299299ZUTT = 14.4
    300300ZUTQ = 10.0
    301 !
     301
    302302ZCDIRU = ZCOEFU(1) + 2.0*ZCOEFU(2)*ZUTU + 3.0*ZCOEFU(3)*ZUTU**2   &
    303303                   + 4.0*ZCOEFU(4)*ZUTU**3 + 5.0*ZCOEFU(5)*ZUTU**4
    … …  
    305305                   + 4.0*ZCOEFT(4)*ZUTT**3
    306306ZCDIRQ = ZCOEFQ(1) + 2.0*ZCOEFQ(2)*ZUTQ
    307 !
     307
    308308ZORDOU = ZCOEFU(0) + ZCOEFU(1)*ZUTU + ZCOEFU(2)*ZUTU**2 + ZCOEFU(3)*ZUTU**3   &
    309309                   + ZCOEFU(4)*ZUTU**4 + ZCOEFU(5)*ZUTU**5
    … …  
    311311                   + ZCOEFT(4)*ZUTT**4
    312312ZORDOQ = ZCOEFQ(0) + ZCOEFQ(1)*ZUTQ + ZCOEFQ(2)*ZUTQ**2
    313 !
    314 !-------------------------------------------------------------------------------
    315 !
     313
     314!-------------------------------------------------------------------------------
     315
    316316!       1.   AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES.
    317317!       --------------------------------------------------------------------
    318 !
     318
    319319ZDIRCOSZW(:) = 1.0
    320 !
     320
    321321ZETV    = XRV/XRD-1.0   !~0.61 (cf Liu et al. 1979)
    322322ZRDSRV  = XRD/XRV       !~0.622
    … …  
    326326ZBTA    = 16.0
    327327ZGMA    = 7.0           !initially =4.7, modified to 7.0 following G. Caniaux
    328 !
     328
    329329ZP00    = 1013.25E+02
    330 !
     330
    331331PCD  = XUNDEF
    332332PCH  = XUNDEF
    … …  
    339339ZHF  = XUNDEF
    340340ZEF  = XUNDEF
    341 !
     341
    342342PSFTH  = XUNDEF
    343343PSFTQ  = XUNDEF
    … …  
    345345PRESA  = XUNDEF
    346346PRI    = XUNDEF
    347 !
     347
    348348ZTAUR   = 0.0
    349349ZRF     = 0.0
    350350ZEFWEBB = 0.0
    351 !
    352 !-------------------------------------------------------------------------------
    353 !
     351
     352!-------------------------------------------------------------------------------
     353
    354354!       2.   INITIALISATIONS BEFORE ITERATIVE LOOP.
    355355!       -------------------------------------------
    356 !
     356
    357357!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))    !set a minimum value to wind
    358358ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
    … …  
    360360write(*,*) "ZVMOD ",SIZE(ZVMOD)
    361361
    362 !
    363362!       2.0. Radiative fluxes - For warm layer & cool skin
    364 !
     363
    365364!       2.0b. Warm Layer correction
    366 !
     365
    367366!       2.1. Specific humidity at saturation
    368 !
     367
    369368WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
    370369  PQSATA = QSAT_SEAWATER2 (PSST(:),PPS(:),PSSS(:))    !at sea surface
    … …  
    377376!### OLIVIER POUR PRESSION SATURANTE #####
    378377!-------------------------------------------------------------------------------
    379 !
     378
    380379ZFOES  = 1 !PSAT(PT(:))
    381380ZFOES  = 0.98*ZFOES
    382 !
     381
    383382ZWORK1    = ZFOES/PPS
    384383ZWORK2    = XRD/XRV
    … …  
    398397write(*,*) "PQSATA : ",PQSATA
    399398
    400 
    401 !
    402399!*       2.    COMPUTE SATURATION HUMIDITY
    403400!              ---------------------------
    404 !
     401
    405402!PQSAT  = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1)
    406403!ZQSATA = ZWORK2A*ZWORK1A / (1.+(ZWORK2A-1.)*ZWORK1A)
    407404ZQSATA = QSAT_SEAWATER (PTA(:),PPA(:))            !at sea surface
    408405
    409 
    410 !
    411406!       2.2. Gradients at the air-sea interface
    412 !
     407
    413408ZDU(:) = ZVMOD(:)               !one assumes u is measured / sea surface current
    414409ZDT(:) = PTA(:)/PEXNA(:)-PSST(:)/PEXNS(:)
    … …  
    418413write(*,*) "PQSAT",PQSAT(:)
    419414write(*,*) "ZDQ",ZDQ(:)
    420 !
     415
    421416!       2.3. Latent heat of vaporisation
    422 !
     417
    423418ZLVA(:) = XLVTT+(XCPV-XCL)*(PTA (:)-XTT)                !of pure water at atm level
    424419ZLVS(:) = XLVTT+(XCPV-XCL)*(PSST(:)-XTT)                !of pure water at sea surface
    … …  
    433428  ZLVS(:) = ZLVS(:)*(1.0-1.00472E-3*PSSS(:))            !of seawater at sea surface
    434429ENDWHERE
    435 !
     430
    436431!       2.4. Specific heat of moist air (Businger 1982)
    437 !
     432
    438433!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
    439434ZCPA(:) = XCPD
    440 !
     435
    441436!       2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
    442 !
     437
    443438ZVISA(:) = 1.326E-05*(1.0+6.542E-03*(PTA(:)-XTT)+8.301E-06*(PTA(:)-XTT)**2   &
    444439           -4.84E-09*(PTA(:)-XTT)**3)
    445 !
     440
    446441!       2.4c Coefficients for warm layer and/or cool skin correction
    447 !
     442
    448443!       2.5. Initial guess
    449 !
     444
    450445ZDDU(:) = ZDU(:)
    451446ZDDT(:) = ZDT(:)
    … …  
    459454write(*,*) "ZDDT ",ZDDT
    460455
    461 !
    462456JCV (:) = -1
    463457ZUSR(:) = 0.04*ZDDU(:)
    … …  
    468462ZDELTAQ10N(:) = ZDDQ(:)
    469463JITER(:) = 99
    470 !
     464
    471465! In the following, we suppose that Richardson number PRI < XRIMAX
    472466! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
    473467!-------------------------------------------------------------------------------
    474 !
     468
    475469!       3.   ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
    476470!       ------------------------------------------
    477 !
     471
    478472DO JJ=1,NITERFL
    479473  DO JLON=1,SIZE(PTA)
    480 !
     474
    481475  IF (JCV(JLON) == -1) THEN
    482476    ZUSR0(JLON)=ZUSR(JLON)
    … …  
    492486    ZDTSTO(JLON) = ZDELTAT10N(JLON)
    493487    ZDQSTO(JLON) = ZDELTAQ10N(JLON)
    494 !
     488
    495489!       3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
    496 !
     490
    497491    IF (ZDELTAU10N(JLON) <= ZUTU) THEN
    498492      ZPARUN(JLON) = ZCOEFU(0) + ZCOEFU(1)*ZDELTAU10N(JLON)      &
    … …  
    505499    ENDIF
    506500    PCDN(JLON) = (ZPARUN(JLON)/ZDELTAU10N(JLON))**2
    507 !
     501
    508502!       3.2. Neutral parameter for temperature (ECUME_V6 formulation)
    509 !
     503
    510504    IF (ZDELTAU10N(JLON) <= ZUTT) THEN
    511505      ZPARTN(JLON) = ZCOEFT(0) + ZCOEFT(1)*ZDELTAU10N(JLON)      &
    … …  
    516510      ZPARTN(JLON) = ZCDIRT*(ZDELTAU10N(JLON)-ZUTT) + ZORDOT
    517511    ENDIF
    518 !
     512
    519513!       3.3. Neutral parameter for humidity (ECUME_V6 formulation)
    520 !
     514
    521515    IF (ZDELTAU10N(JLON) <= ZUTQ) THEN
    522516      ZPARQN(JLON) = ZCOEFQ(0) + ZCOEFQ(1)*ZDELTAU10N(JLON)      &
    … …  
    525519      ZPARQN(JLON) = ZCDIRQ*(ZDELTAU10N(JLON)-ZUTQ) + ZORDOQ
    526520    ENDIF
    527 !
     521
    528522!       3.4. Scaling parameters U*, T*, Q*
    529 !
     523
    530524    ZUSR(JLON) = ZPARUN(JLON)
    531525    ZTSR(JLON) = ZPARTN(JLON)*ZDELTAT10N(JLON)/ZDELTAU10N(JLON)
    532526    ZQSR(JLON) = ZPARQN(JLON)*ZDELTAQ10N(JLON)/ZDELTAU10N(JLON)
    533 !
     527
    534528!       3.4b Gustiness factor (Deardorff 1970)
    535 !
     529
    536530!       3.4c Cool skin correction
    537 !
     531
    538532!       3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
    539 !
     533
    540534! For U
    541535    ZLMOU = PUREF(JLON)*XG*XKARMAN*(ZTSR(JLON)/PTA(JLON)   &
    … …  
    545539    ZLMOU = MAX(MIN(ZLMOU,ZLMOMAX),ZLMOMIN)
    546540    ZLMOT = MAX(MIN(ZLMOT,ZLMOMAX),ZLMOMIN)
    547 !
     541
    548542!       3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
    549543!            Modified to include convective form following Fairall (unpublished)
    550 !
     544
    551545! For U
    552546    IF (ZLMOU == 0.0) THEN
    … …  
    581575    ENDIF
    582576    ZPSIT(JLON) = ZPSI_T
    583 !
     577
    584578!       3.7. Update air-sea gradients
    585 !
     579
    586580    ZDDU(JLON) = ZDU(JLON)
    587581    ZDDT(JLON) = ZDT(JLON)
    … …  
    601595    ZDELTAQ10N(JLON) = SIGN(MAX(ABS(ZDELTAQ10N(JLON)),10.0*ZDQSR0),   &
    602596                            ZDELTAQ10N(JLON))
    603 !
     597
    604598!       3.8. Test convergence for U*, T*, Q*
    605 !
     599
    606600    IF (ABS(ZUSR(JLON)-ZUSR0(JLON)) < ZDUSR0 .AND.   &
    607601        ABS(ZTSR(JLON)-ZTSR0(JLON)) < ZDTSR0 .AND.   &
    … …  
    612606    JITER(JLON) = JJ
    613607  ENDIF
    614 !
     608
    615609  ENDDO
    616610ENDDO
    617 !
    618 !-------------------------------------------------------------------------------
    619 !
     611
     612!-------------------------------------------------------------------------------
     613
    620614!       4.   COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
    621615!       ---------------------------------------------------------------
    622 !
     616
    623617DO JLON=1,SIZE(PTA)
    624 !
     618
    625619!       4.1. Surface turbulent fluxes
    626620!            (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
    627 !
     621
    628622  ZTAU(JLON) = -PRHOA(JLON)*ZUSR(JLON)**2
    629623  ZHF(JLON)  = -PRHOA(JLON)*ZCPA(JLON)*ZUSR(JLON)*ZTSR(JLON)
    … …  
    634628  write(*,*) "LAT = ",ZEF(JLON)
    635629
    636 !
    637630!       4.2. Exchange coefficients PCD, PCH, PCE
    638 !
     631
    639632  PCD(JLON) = (ZUSR(JLON)/ZDDU(JLON))**2
    640633  PCH(JLON) = (ZUSR(JLON)*ZTSR(JLON))/(ZDDU(JLON)*ZDDT(JLON))
    … …  
    645638  write(*,*) "ZQSR = ",ZQSR(JLON)
    646639
    647 !
    648640!       4.3. Stochastic perturbation of turbulent fluxes
    649 !
     641
    650642!  IF( OPERTFLUX )THEN
    651643!    ZTAU(JLON) = ZTAU(JLON)* ( 1. + PPERTFLUX(JLON) / 2. )
    … …  
    653645!    ZEF (JLON) = ZEF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
    654646!  ENDIF
    655 !
     647
    656648ENDDO
    657 !
    658 !-------------------------------------------------------------------------------
    659 !
     649
     650!-------------------------------------------------------------------------------
     651
    660652!       5.   COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
    661653!            (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
    662654!       -----------------------------------------------------
    663 !
     655
    664656IF (OPRECIP) THEN
    665657  DO JLON=1,SIZE(PTA)
    666 !
     658
    667659!       5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
    668 !
     660
    669661! See pp3752 in FBR96.
    670662    ZTAUR(JLON) = -0.85*PRAIN(JLON)*PVMOD(JLON)
    671 !
     663
    672664!       5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
    673 !
     665
    674666! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
    675667! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity
    676668! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
    677669! as wet-bulb factor (Eq.11 in GoF95).
    678 !
     670
    679671    ZTAC   = PTA(JLON)-XTT
    680672    ZCPWA  = 4217.51 -3.65566*ZTAC +0.1381*ZTAC**2       &
    … …  
    688680    ZRF(JLON) = PRAIN(JLON)*ZCPWA*ZBULB*((PSST(JLON)-PTA(JLON))   &
    689681                +(PQSATA(JLON)-PQA(JLON))*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
    690 !
     682
    691683  ENDDO
    692684ENDIF
    693 !
    694 !-------------------------------------------------------------------------------
    695 !
     685
     686!-------------------------------------------------------------------------------
     687
    696688!       6.   COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
    697689!       ------------------------------------------------------------------------
    698 !
     690
    699691! See Eq.21 and Eq.22 in FBR96.
    700692IF (OPWEBB) THEN
    … …  
    705697  ENDDO
    706698ENDIF
    707 !
    708 !-------------------------------------------------------------------------------
    709 !
     699
     700!-------------------------------------------------------------------------------
     701
    710702!       7.   FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS.
    711703!       ---------------------------------------------------------------
    712 !
     704
    713705!       7.1. Richardson number
    714 !
     706
    715707! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,PQA,   &
    716708!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI)
    717 !
     709
    718710!       7.2. Friction velocity which contains correction due to rain
    719 !
     711
    720712ZUSTAR2(:) = -(ZTAU(:)+ZTAUR(:))/PRHOA(:)       !>>0 as ZTAU<<0 & ZTAUR<=0
    721 !
     713
    722714IF (OPRECIP) THEN
    723715  PCD(:) = ZUSTAR2(:)/ZDDU(:)**2
    724716ENDIF
    725 !
     717
    726718PUSTAR(:) = SQRT(ZUSTAR2(:))                    !>>0
    727 !
     719
    728720!       7.3. Aerodynamical conductance and resistance
    729 !
     721
    730722ZAC  (:) = PCH(:)*ZDDU(:)
    731723PRESA(:) = 1.0/ZAC(:)
    732 !
     724
    733725!       7.4. Total surface fluxes
    734 !
     726
    735727PSFTH(:) =  ZHF(:)+ZRF(:)
    736728PSFTQ(:) = (ZEF(:)+ZEFWEBB(:))/ZLVS(:)
    737 !
     729
    738730!       7.5. Charnock number
    739 !
     731
    740732IF (CCHARNOCK == 'OLD') THEN
    741733  ZCHARN(:) = XVCHRNK
    … …  
    743735  ZCHARN(:) = MIN(0.018,MAX(0.011,0.011+(0.007/8.0)*(ZDDU(:)-10.0)))
    744736ENDIF
    745 !
     737
    746738!       7.6. Roughness lengths Z0 and Z0H over sea
    747 !
     739
    748740!IF (KZ0 == 0) THEN      ! ARPEGE formulation
    749741!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + XVZ0CM*PCD(:)/PCDN(:)
    … …  
    770762       PCE,&
    771763       PCH]
    772 !
    773 !
     764
     765
    774766!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',1,ZHOOK_HANDLE)
    775 !
     767
    776768!-------------------------------------------------------------------------------
    777769   END SUBROUTINE ECUMEV6_FLUX
Note: See TracChangeset for help on using the changeset viewer.