Changeset 6176


Ignore:
Timestamp:
Apr 18, 2026, 10:17:55 PM (2 weeks ago)
Author:
evignon
Message:

modification de lmdz_cloud_optics pour inclure une option deboguee
du calcul du rayon optique des gouttelettes en presence d'aerosols.
On en profite pour retructurer, nettoyer et commenter la routine.
Convergence assuree en debug et prod sans aerosols.
Convergence assuree en debug avec aerosols, mais perdue en prod si
aerosols et ok_cdnc = y en raison de refactorisations

Location:
LMDZ6/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_cloud_optics_prop.f90

    r6154 r6176  
    1010CONTAINS
    1111
     12!=================================================================================================       
    1213SUBROUTINE cloud_optics_prop_post()
    1314  USE lmdz_cloud_optics_prop_ini, ONLY: novlp
    … …  
    3435END SUBROUTINE cloud_optics_prop_post
    3536
     37
     38!=================================================================================================
    3639SUBROUTINE cloud_optics_prop(klon, klev, paprs, pplay, temp, radocond, picefra, pclc, &
    37     pcltau, pclemi, pch, pcl, pcm, pct, radocondwp, xflwp, xfiwp, xflwc, xfiwc, &
     40    pcldtau, pclemi, pch, pcl, pcm, pct, radocondwp, xflwp, xfiwp, xflwc, xfiwc, &
    3841    mass_solu_aero, mass_solu_aero_pi, pcldtaupi, distcltop, temp_cltop, re, fl, reliq, reice, &
    3942    reliq_pi, reice_pi, scdnc, cdnc, cdnc_pi, cldncl, reffclwtop, lcc, reffclws, &
    … …  
    5962  USE lmdz_cloud_optics_prop_ini , ONLY : rei_coef, rei_min_temp
    6063  USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
    61   USE lmdz_cloud_optics_prop_ini , ONLY : first
     64  USE lmdz_cloud_optics_prop_ini , ONLY : first, iflag_cdreff
    6265 
    6366
    … …  
    7780  !          "atelier optics of clouds": careful reading and comments
    7881  !                                   things to address -> see ?aoc
     82  !          E. Vignon: comments + cleaning + structuration
    7983  !
    8084  ! Aim: compute condensates' optical properties
    81   !      (cloud optical depth,) and emissivity ?aoc
     85  !      (cloud optical depth,) and emissivity
    8286  ! ======================================================================
    8387 
    … …  
    126130 
    127131  ! diagnostics or properties if one uses oldrad
    128   REAL, INTENT(OUT) :: pcltau(klon, klev) ! cloud optical depth [m]
     132  REAL, INTENT(OUT) :: pcldtau(klon, klev) ! cloud optical depth [m]
    129133  REAL, INTENT(OUT) :: pclemi(klon, klev) ! cloud emissivity [-]
    130   REAL, INTENT(OUT) :: pcldtaupi(klon, klev) ! pre-industrial value of cloud optical thickness, ie.
    131                                              ! values of optical thickness that does not account
     134  REAL, INTENT(OUT) :: pcldtaupi(klon, klev) ! pre-industrial value of cloud optical depth, ie.
     135                                             ! values of optical depth that does not account
    132136                                             ! for aerosol effects on cloud droplet radius [m]
    133137  REAL, INTENT(OUT) :: scdnc(klon, klev)   ! cloud droplet number concentration, mean over the whole mesh [m-3]
    … …  
    158162  INTEGER flag_max
    159163
    160   ! threshold PARAMETERs
     164  ! threshold parameters
    161165  REAL phase3d(klon, klev)
    162166  REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
    … …  
    168172  REAL rel, tc, rei, iwc, dei, deimin, deimax
    169173  REAL k_ice
     174  REAL rhol ! density of liquid water in kg/m3
    170175  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
    171176
    172   ! IM cf. CR:parametres supplementaires
    173177  REAL dzfice(klon,klev)
    174178  REAL zclear(klon)
    … …  
    177181  REAL zcloudm(klon)
    178182  REAL zcloudl(klon)
    179   REAL rhodz(klon, klev) !--rho*dz pour la couche
    180   REAL zrho(klon, klev) !--rho pour la couche
    181   REAL dh(klon, klev) !--dz pour la couche
    182   REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
    183   REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
     183  REAL rhodz(klon, klev) !--rho*dz of layer
     184  REAL zrho(klon, klev) !--rho of layer
     185  REAL dh(klon, klev) !--dz of layer
     186  REAL rad_chaud(klon, klev) !--radius for warm ("chaud") liquid clouds
     187  REAL rad_chaud_pi(klon, klev) !--radius for warm ("chaud") liquid clouds, pre-industrial
    184188  REAL zflwp_var, zfiwp_var
    185189  REAL d_rei_dt
    186 
    187 
    188   ! FH : 2011/05/24
    189   ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
    190   ! to be used for a temperature in celcius T(°C) < 0
    191   ! rei=rei_min for T(°C) < -81.4
    192   ! Calcul de la pente de la relation entre rayon effective des cristaux
    193   ! et la température Pour retrouver les résultats numériques de la version d'origine,
    194   ! on impose 0.71 quand on est proche de 0.71
    195   d_rei_dt = (rei_max-rei_min)/81.4
    196   IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
    197 
    198   ! Calculer l'epaisseur optique et l'emmissivite des nuages
    199   ! IM inversion des DO
    200 
    201   xflwp = 0.D0
    202   xfiwp = 0.D0
    203   xflwc = 0.D0
    204   xfiwc = 0.D0
    205 
    206   reliq = 0.
    207   reice = 0.
    208   reliq_pi = 0.
    209   reice_pi = 0.
    210 
     190!=====================================================================================
     191
     192! Initialisation
     193!===============
     194  xflwp(:) = 0.D0
     195  xfiwp(:) = 0.D0
     196  xflwc(:,:) = 0.D0
     197  xfiwc(:,:) = 0.D0
     198  radocondwp(:) = 0.
     199  rhol=1000.0
     200
     201  reliq(:,:) = 0.
     202  reice(:,:) = 0.
     203  reliq_pi(:,:) = 0.
     204  reice_pi(:,:) = 0.
     205
     206! Preliminary calculations
     207!==========================
     208
     209
     210! for 'old' (pre-CMIP7) large-scale condensation scheme
     211! (firstilp), ice fraction is recomputed here
    211212  IF ((.NOT. ok_new_lscp) .AND. iflag_t_glace.EQ.0) THEN
    212213    DO k = 1, klev
    213214      DO i = 1, klon
    214         ! -layer calculation
    215         rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
    216         zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
    217         dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
    218         ! -Fraction of ice in cloud using a linear transition
    219215        icefrac_optics(i, k) = 1.0 - (temp(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
    220216        icefrac_optics(i, k) = min(max(icefrac_optics(i,k),0.0), 1.0)
    221         ! -IM Total Liquid/Ice water content
    222         xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
    223         xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
    224       ENDDO
    225     ENDDO
    226   ELSE ! of IF (iflag_t_glace.EQ.0)
    227     DO k = 1, klev
    228 
    229 
    230       DO i = 1, klon
    231 
     217      ENDDO
     218    ENDDO
     219
     220  ELSE
     221! ice fraction directly from the lscp param, except for convective grid point           
     222! where we take the temperature-dependent icefrac_optics computed in lmdz_call_cloud_optics
     223
     224   DO k = 1, klev
     225      DO i = 1, klon
    232226        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
    233         ! EV: take the ice fraction directly from the lscp code
    234         ! consistent only for non convective grid points
    235         ! critical for mixed phase clouds
    236             icefrac_optics(i,k)=picefra(i,k)
     227           icefrac_optics(i,k)=picefra(i,k) ! picefra is the ice fraction computed in lscp
    237228        ENDIF
    238 
    239         ! -layer calculation
    240         rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
    241         zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
    242         dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
    243         ! -IM Total Liquid/Ice water content
    244         xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
    245         xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
    246229      ENDDO
    247230    ENDDO
    248231  ENDIF
    249232
    250 
    251 
    252 
    253 
     233! computation of layers'mass and water contents for radiation
     234DO k = 1,klev
     235  DO i = 1,klon
     236    rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
     237    zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
     238    dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
     239    ! -Liquid/Ice water specific content, mesh averaged
     240    xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
     241    xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
     242   ! vertically integrated contents (water paths)
     243    xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
     244    xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
     245    radocondwp(i) = radocondwp(i) + radocond(i, k)*rhodz(i, k)
     246  END DO
     247END DO
     248
     249
     250
     251! Computation of cloud droplet effective radius
     252!===============================================
     253! There 2 options to compute cloud droplet effective radius:
     254! if ok_cdnc: we compute cloud radius as a function of an in-cloud number concentration of cloud droplets (cdnc)
     255! depending on the concentration of soluble aerosols
     256! otherwise, cloud droplet radius is fixed to a constant value (in fact one constant for the first three layers, another
     257! one for layers above
     258
     259
     260! if ok_cdnc, let's first compute the in-cloud droplet number concentration
    254261
    255262  IF (ok_cdnc) THEN
    256 
    257     ! --we compute cloud properties as a function of the aerosol load
    258 
     263   
     264    ! Pre-industrial cloud droplet concentrations
    259265    DO k = 1, klev
    260266      DO i = 1, klon
    … …  
    262268        ! Cloud droplet number concentration (CDNC) is restricted
    263269        ! to be within [20, 1000 cm^3]
    264 
    265         ! --pre-industrial case
    266270        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
    267271          1.E-4))/log(10.))*1.E6 !-m-3
    … …  
    271275    ENDDO
    272276
    273     !--flag_aerosol=7 => MACv2SP climatology
    274     !--in this case there is an enhancement factor
    275     IF (flag_aerosol .EQ. 7) THEN
    276 
    277       !--present-day
     277   ! Present day cloud droplet concentrations
     278   IF (flag_aerosol .EQ. 7) THEN
     279
     280    ! flag_aerosol=7 => MACv2SP climatology
     281    ! in this case we apply an enhancement factor dNovrN on the pi value
    278282      DO k = 1, klev
    279283        DO i = 1, klon
    … …  
    282286      ENDDO
    283287
    284     !--standard case
    285     ELSE
    286 
     288   ELSE
     289   ! standard configuration using the same formula as above but
     290   ! considering present day aerosols concentrations
    287291      DO k = 1, klev
    288292        DO i = 1, klon
    289 
    290           ! Formula "D" of Boucher and Lohmann, Tellus, 1995
    291           ! Cloud droplet number concentration (CDNC) is restricted
    292           ! to be within [20, 1000 cm^3]
    293 
    294           ! --present-day case
    295293          cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
    296294            1.E-4))/log(10.))*1.E6 !-m-3
    297295          cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
    298 
    299296        ENDDO
    300297      ENDDO
    … …  
    302299    ENDIF !--flag_aerosol
    303300
    304     !--computing cloud droplet size
     301 ENDIF
     302
     303! Cloud droplet effective radius calculation
     304! controled by ok_cdnc and iflag_cdreff
     305
     306 IF (ok_cdnc .AND. iflag_cdreff .EQ. 0) THEN
     307    ! aerosol-aware formulation used until CMIP6
     308    ! Warning! the formula is erroneous as it
     309    ! considers the liq+ice water mass and
     310    ! the mesh-averaged (and not the in-cloud)
     311    ! contents
    305312    DO k = 1, klev
    306313      DO i = 1, klon
     314
     315        ! --pre-industrial case
     316        rad_chaud_pi(i, k) = 1.1*((radocond(i,k)*pplay(i, &
     317          k)/(rd*temp(i,k)))/(4./3.*rpi*rhol*cdnc_pi(i,k)))**(1./3.)
     318        rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
    307319
    308320        ! --present-day case
    309321        rad_chaud(i, k) = 1.1*((radocond(i,k)*pplay(i, &
    310           k)/(rd*temp(i,k)))/(4./3*rpi*1000.*cdnc(i,k)))**(1./3.)
     322          k)/(rd*temp(i,k)))/(4./3*rpi*rhol*cdnc(i,k)))**(1./3.)
    311323        rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
    312324
     325        END DO
     326    ENDDO
     327
     328 ELSE IF (ok_cdnc .AND. iflag_cdreff .EQ. 1) THEN
     329     ! similar as above but debugged
     330
     331    DO k = 1, klev
     332      DO i = 1, klon
     333
    313334        ! --pre-industrial case
    314         rad_chaud_pi(i, k) = 1.1*((radocond(i,k)*pplay(i, &
    315           k)/(rd*temp(i,k)))/(4./3.*rpi*1000.*cdnc_pi(i,k)))**(1./3.)
    316         rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
    317 
    318         ! --pre-industrial case
    319         ! --liquid/ice cloud water paths:
    320         IF (pclc(i,k)<=seuil_neb) THEN
    321 
    322           pcldtaupi(i, k) = 0.0
    323 
    324         ELSE
    325 
    326           zflwp_var = 1000.*(1.-icefrac_optics(i,k))*radocond(i, k)/pclc(i, k)* &
    327             rhodz(i, k)
    328           zfiwp_var = 1000.*icefrac_optics(i, k)*radocond(i, k)/pclc(i, k)*rhodz(i, k)
    329           ! Calculation of ice cloud effective radius in micron
    330           IF (iflag_rei .EQ. 2) THEN
    331             ! in-cloud ice water content in g/m3
    332             iwc = icefrac_optics(i,k) * radocond(i,k) / pclc(i,k) * zrho(i,k) * 1000.
    333             ! this formula is a simplified version of the Sun 2001 one (as in the IFS model,
    334             ! and which is activated for iflag_rei = 1).
    335             ! In particular, the term in temperature**2 has been simplified.
    336             ! The new coefs are tunable, and are by default set so that the results fit well
    337             ! to the Sun 2001 formula
    338             ! By default, rei_coef = 2.4 and rei_min_temp = 175.
    339             ! The ranges of these parameters are determined so that the RMSE between this
    340             ! formula and the one from Sun 2001 is less than 4 times the minimum RMSE
    341             ! The ranges are [1.9, 2.9] for rei_coef and [160., 185.] for rei_min_temp
    342             dei = rei_coef * (iwc**0.2445) * ( temp(i,k) - rei_min_temp )
    343             ! we clip the results
    344             deimin = 20.
    345             deimax = 155.
    346             dei = MIN(MAX(dei, deimin), deimax)
    347             ! formula to convert effective diameter in effective radius
    348             rei = 3. * SQRT(3.) / 8. * dei
    349           ELSEIF (iflag_rei .EQ. 1) THEN
    350             ! when we account for precipitation in the radiation scheme,
    351             ! It is recommended to use the rei formula from Sun and Rikkus 1999 with a revision
    352             ! from Sun 2001 (as in the IFS model)
    353             iwc=icefrac_optics(i, k)*radocond(i, k)/pclc(i,k)*zrho(i,k)*1000. !in cloud ice water content in g/m3
    354             dei=(1.2351+0.0105*(temp(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
    355                & 0.7957*(iwc**0.2535)*(temp(i,k)-83.15))
    356             !deimax=155.0
    357             !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
    358             !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
    359             deimax=rei_max*2.0
    360             deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
    361             dei=min(dei,deimax)
    362             dei=max(dei,deimin)
    363             rei=3.*sqrt(3.)/8.*dei
    364            ELSE
    365             ! Default
    366             ! for ice clouds: as a function of the ambiant temperature
    367             ! [formula used by Iacobellis and Somerville (2000), with an
    368             ! asymptotical value of 3.5 microns at T<-81.4 C added to be
    369             ! consistent with observations of Heymsfield et al. 1986]:
    370             ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
    371             ! rei_max=61.29
    372             tc = temp(i, k) - 273.15
    373             rei = d_rei_dt*tc + rei_max
    374             IF (tc<=-81.4) rei = rei_min
    375            ENDIF
    376 
    377           ! -- cloud optical thickness :
    378           ! [for liquid clouds, traditional formula,
    379           ! for ice clouds, Ebert & Curry (1992)]
    380 
    381           IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    382           pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
    383             zfiwp_var*(3.448E-03+2.431/rei)
    384 
    385         ENDIF
    386 
    387       ENDDO
    388     ENDDO
    389 
    390   ELSE !--not ok_cdnc
    391 
    392     ! -prescribed cloud droplet radius
    393 
     335        rad_chaud_pi(i, k) = (xflwc(i,k)/max(pclc(i,k),seuil_neb)*pplay(i, &
     336          k)/(rd*temp(i,k))/(4./3.*rpi*rhol*cdnc_pi(i,k)))**(1./3.)
     337        rad_chaud_pi(i, k) = min(max(rad_chaud_pi(i,k)*1.E6, 5.),100.)
     338
     339        ! --present-day case
     340        rad_chaud(i, k) = (xflwc(i,k)/max(pclc(i,k),seuil_neb)*pplay(i, &
     341          k)/(rd*temp(i,k))/(4./3*rpi*rhol*cdnc(i,k)))**(1./3.)
     342        rad_chaud(i, k) = min(max(rad_chaud(i,k)*1.E6, 5.),100.)
     343
     344        END DO
     345    ENDDO
     346
     347 ELSE ! ok_cdnc
     348
     349    ! -fixed (prescribed) cloud droplet effective radius values
    394350    DO k = 1, min(3, klev)
    395351      DO i = 1, klon
    … …  
    404360      ENDDO
    405361    ENDDO
    406 
    407   ENDIF !--ok_cdnc
    408 
    409   ! --computation of cloud optical depth and emissivity
    410   ! --in the general case
     362 
     363 ENDIF ! ok_cdnc
     364
     365  ! For output diagnostics of cloud droplet effective radius [um]
     366  ! we multiply here with f * xl (fraction of liquid water
     367  ! clouds in the grid cell) to avoid problems in the averaging of the output.
     368  ! The actual cloud droplet effective radius can be computed from outputs as re/fl
    411369
    412370  DO k = 1, klev
    413371    DO i = 1, klon
    414372
     373      reliq(i, k) = rad_chaud(i, k)
     374      reliq_pi(i, k) = rad_chaud_pi(i, k)     
     375
    415376      IF (pclc(i,k)<=seuil_neb) THEN
    416 
    417         ! effective cloud droplet radius (microns) for liquid water clouds:
    418         ! For output diagnostics cloud droplet effective radius [um]
    419         ! we multiply here with f * xl (fraction of liquid water
    420         ! clouds in the grid cell) to avoid problems in the averaging of the
    421         ! output.
    422         ! In the output of IOIPSL, derive the REAL cloud droplet
    423         ! effective radius as re/fl
    424377
    425378        fl(i, k) = seuil_neb*(1.-icefrac_optics(i,k))
    426379        re(i, k) = rad_chaud(i, k)*fl(i, k)
    427         rel = 0.
    428         rei = 0.
    429         pclc(i, k) = 0.0
    430         pcltau(i, k) = 0.0
    431         pclemi(i, k) = 0.0
    432380
    433381      ELSE
    434 
    435         ! -- liquid/ice cloud water paths:
    436 
    437         zflwp_var = 1000.*(1.-icefrac_optics(i,k))*radocond(i, k)/pclc(i, k)*rhodz(i, k)
    438         zfiwp_var = 1000.*icefrac_optics(i, k)*radocond(i, k)/pclc(i, k)*rhodz(i, k)
    439 
    440         ! effective cloud droplet radius (microns) for liquid water clouds:
    441         ! For output diagnostics cloud droplet effective radius [um]
    442         ! we multiply here with f Effective radius of cloud droplet at top of cloud (m)* xl (fraction of liquid water
    443         ! clouds in the grid cell) to avoid problems in the averaging of the
    444         ! output.
    445         ! In the output of IOIPSL, derive the REAL cloud droplet
    446         ! effective radius as re/fl
    447382
    448383        fl(i, k) = pclc(i, k)*(1.-icefrac_optics(i,k))
    449384        re(i, k) = rad_chaud(i, k)*fl(i, k)
    450 
    451         rel = rad_chaud(i, k)
     385 
     386       ENDIF
     387
     388    END DO
     389  END DO
     390
     391! Computation of ice crystal effective radius
     392!=============================================
     393
     394  DO k = 1, klev
     395    DO i = 1, klon
     396
     397      IF (pclc(i,k)>seuil_neb) THEN
    452398
    453399        ! Calculation of ice cloud effective radius in micron
    454 
    455400
    456401        IF (iflag_rei .EQ. 2) THEN
    … …  
    498443            ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
    499444            ! rei_max=61.29
     445           
     446            ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
     447            ! to be used for a temperature in celcius T(°C) < 0
     448            ! rei=rei_min for T(°C) < -81.4
     449            ! Computation of slope for ice crystal effective radius=f(T)
     450            ! For consistency with original version, we impoe 0.71 when
     451            ! close to this value
     452            d_rei_dt = (rei_max-rei_min)/81.4
     453            IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71         
     454           
    500455            tc = temp(i, k) - 273.15
    501456            rei = d_rei_dt*tc + rei_max
    502457            IF (tc<=-81.4) rei = rei_min
     458
    503459        ENDIF
     460
     461      ELSE ! pclc > rneb
     462
     463           rei = 0.
     464
     465      ENDIF
     466
     467      reice(i, k) = rei
     468      ! same radius for pi conditions
     469      reice_pi(i,k) = reice(i,k)
     470
     471    ENDDO
     472  ENDDO
     473
     474 
     475
     476! Computation of cloud optical depth and emissivity
     477!===================================================
     478
     479  DO k = 1, klev
     480    DO i = 1, klon
     481
     482      IF (pclc(i,k)<=seuil_neb) THEN
     483
     484        rel = 0.
     485        rei = 0.
     486        pcldtau(i, k) = 0.0
     487        pcldtaupi(i, k) = 0.0
     488        pclemi(i, k) = 0.0
     489
     490      ELSE
     491
     492        rel = rad_chaud(i, k)
     493        rei = reice(i,k)
     494        ! mass of in-cloud condensates in g
     495        zflwp_var = 1000.*(1.-icefrac_optics(i,k))*radocond(i, k)/pclc(i, k)*rhodz(i, k)
     496        zfiwp_var = 1000.*icefrac_optics(i, k)*radocond(i, k)/pclc(i, k)*rhodz(i, k)
     497        IF (zflwp_var==0.) rel = 1.
     498        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    504499        ! -- cloud optical thickness :
    505500        ! [for liquid clouds, traditional formula,
    506         ! for ice clouds, Ebert & Curry (1992)]
    507 
    508         IF (zflwp_var==0.) rel = 1.
    509         IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    510         pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
    511           rei)
    512 
     501        ! for ice clouds, Ebert & Curry (1992)]       
     502        pcldtau(i, k) = 3.0/2.0*(zflwp_var/rel) &
     503                    + zfiwp_var*(3.448E-03+2.431/ rei)
     504        IF (ok_cdnc) THEN
     505            pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
     506                            zfiwp_var*(3.448E-03+2.431/rei)
     507        ELSE
     508        ! -- if cloud droplet radius is fixed, plcdtaupi = plcdtau
     509            pcldtaupi(i, k) = pcldtau(i,k)
     510        END IF
    513511        ! -- cloud infrared emissivity:
    514         ! [the broadband infrared absorption coefficient is PARAMETERized
    515         ! as a function of the effective cld droplet radius]
     512        ! [the broadband infrared absorption coefficient is parameterized
     513        ! as a function of the effective cloud droplet radius]
    516514        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
    517 
    518515        k_ice = k_ice0 + 1.0/rei
    519 
    520516        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
    521517
    522518      ENDIF
    523519
    524       reice(i, k) = rei
    525 
    526       xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
    527       xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
    528 
    529520    ENDDO
    530521  ENDDO
    531522
    532   ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
    533 
    534   IF (.NOT. ok_cdnc) THEN
    535     DO k = 1, klev
    536       DO i = 1, klon
    537         pcldtaupi(i, k) = pcltau(i, k)
    538         reice_pi(i, k) = reice(i, k)
    539       ENDDO
    540     ENDDO
    541   ENDIF
    542 
    543   DO k = 1, klev
    544     DO i = 1, klon
    545       reliq(i, k) = rad_chaud(i, k)
    546       reliq_pi(i, k) = rad_chaud_pi(i, k)
    547       reice_pi(i, k) = reice(i, k)
    548     ENDDO
    549   ENDDO
    550 
    551   ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
    552   ! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
    553   ! initialisations
     523
     524! Cloud cover calculation
     525! ========================
    554526
    555527  DO i = 1, klon
    … …  
    562534    pcm(i) = 1.0
    563535    pcl(i) = 1.0
    564     radocondwp(i) = 0.0
     536
    565537  ENDDO
    566538
    567   ! --calculation of liquid water path
    568 
    569   DO k = klev, 1, -1
    570     DO i = 1, klon
    571       radocondwp(i) = radocondwp(i) + radocond(i, k)*rhodz(i, k)
    572     ENDDO
    573   ENDDO
    574 
    575   ! --calculation of cloud properties with cloud overlap
    576   ! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
    577   ! !novlp=1: max-random
    578   ! !novlp=2: maximum
    579   ! !novlp=3: random
     539  DO k=1, klev
     540     DO i=1,klon
     541       IF (pclc(i,k) <= seuil_neb) THEN
     542            pclc(i,k) = 0.
     543       END IF
     544     END DO
     545  END DO
     546
     547  ! Cloud overlap approximation
     548  ! choix made through radopt.h
     549  ! novlp=1: max-random
     550  ! novlp=2: maximum
     551  ! novlp=3: random
    580552
    581553
    … …  
    638610  ENDDO
    639611
    640   ! ========================================================
    641   ! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
    642   ! ========================================================
     612  ! Diagnostics computation for CMIP protocol
     613  ! =========================================
    643614  ! change by Nicolas Yan (LSCE)
    644615  ! Cloud Droplet Number Concentration (CDNC) : 3D variable
    … …  
    676647          ! threshold:
    677648
    678         IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
     649        IF (pcldtau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
    679650
    680651          IF (novlp.EQ.2) THEN
    … …  
    720691    ENDDO ! loop over i
    721692
    722     ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
    723     ! REFFCLWS)
     693    ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC, REFFCLWS)
    724694    DO i = 1, klon
    725695      DO k = 1, klev
    … …  
    729699        lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
    730700        lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
    731         !FC pour la glace (CAUSES)
    732701        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*ccwcon(i, k) !  glace convective
    733702        icc3dstra(i, k)= pclc(i, k)*radocond(i, k)*(1-phase3d(i, k))
    734703        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
    735704        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
    736         !FC (CAUSES)
    737705
    738706        ! Compute cloud droplet radius as above in meter
    … …  
    776744        IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
    777745        IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
    778 !FC (CAUSES)
    779746        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
    780747        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
    781 !FC (CAUSES)
    782748      ENDDO
    783749      IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0

  • LMDZ6/trunk/libf/phylmd/lmdz_cloud_optics_prop_ini.f90

    r6140 r6176  
    88  INTEGER, PROTECTED :: iflag_t_glace=0
    99  INTEGER, PROTECTED :: iflag_rei
     10  INTEGER, PROTECTED :: iflag_cdreff
    1011  INTEGER, PROTECTED :: novlp, iflag_ice_thermo
    1112  LOGICAL, PROTECTED :: ok_cdnc
    … …  
    109110    rei_min_temp = 175.
    110111    CALL getin_p('rei_min_temp',rei_min_temp)
    111 
     112    iflag_cdreff = 0
     113    CALL getin_p('iflag_cdreff',iflag_cdreff)
    112114   
    113115  END SUBROUTINE cloud_optics_prop_ini

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

    r6169 r6176  
    10151015            ! ice or snow depending on snow thickness
    10161016            alb_snow = alb_snow - (alb_snow - alb_ice)*EXP(-snow(i)*z1_s)
    1017             ! Effect of clouds (polynomial fit with 83% clouds)
     1017            ! Effect of clouds (polynomial fit with 81% clouds)
    10181018            alb1_new(i) = alb_snow - 0.81*(-0.1010*alb_snow*alb_snow &
    10191019                                          + 0.1933*alb_snow - 0.0148)
Note: See TracChangeset for help on using the changeset viewer.