Ignore:
Timestamp:
Jun 29, 2018, 12:31:11 PM (6 years ago)
Author:
Laurent Fairhead
Message:

First attempt at merging with trunk

Location:
LMDZ6/branches/DYNAMICO-conv
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/newmicro.F90

    r2596 r3356  
    11! $Id$
    22
    3 
    4 
    5 SUBROUTINE newmicro(ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
     3SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
    64    pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
    75    mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
    … …  
    108  USE dimphy
    119  USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
    12     reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra
     10      reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
     11      zfice, dNovrN
    1312  USE phys_state_var_mod, ONLY: rnebcon, clwcon
    1413  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
     14  USE ioipsl_getin_p_mod, ONLY : getin_p
     15  USE print_control_mod, ONLY: lunout
     16
     17
    1518  IMPLICIT NONE
    1619  ! ======================================================================
    … …  
    139142  ! within the grid cell)
    140143
     144  INTEGER flag_aerosol
    141145  LOGICAL ok_cdnc
    142146  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
    … …  
    152156  REAL zrho(klon, klev) !--rho pour la couche
    153157  REAL dh(klon, klev) !--dz pour la couche
    154   REAL zfice(klon, klev)
    155158  REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
    156159  REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
    … …  
    162165  REAL reliq_pi(klon, klev), reice_pi(klon, klev)
    163166
     167  REAL,SAVE :: cdnc_min=-1.
     168  REAL,SAVE :: cdnc_min_m3
     169  !$OMP THREADPRIVATE(cdnc_min,cdnc_min_m3)
     170  REAL,SAVE :: cdnc_max=-1.
     171  REAL,SAVE :: cdnc_max_m3
     172  !$OMP THREADPRIVATE(cdnc_max,cdnc_max_m3)
     173
    164174  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    165175  ! FH : 2011/05/24
    … …  
    173183  ! Pour retrouver les résultats numériques de la version d'origine,
    174184  ! on impose 0.71 quand on est proche de 0.71
     185
     186  if (first) THEN
     187      call getin_p('cdnc_min',cdnc_min)
     188      cdnc_min_m3=cdnc_min*1.E6
     189      IF (cdnc_min_m3<0.) cdnc_min_m3=20.E6 ! astuce pour retrocompatibilite
     190      write(lunout,*)'cdnc_min=', cdnc_min_m3/1.E6
     191      call getin_p('cdnc_max',cdnc_max)
     192      cdnc_max_m3=cdnc_max*1.E6
     193      IF (cdnc_max_m3<0.) cdnc_max_m3=1000.E6 ! astuce pour retrocompatibilite
     194      write(lunout,*)'cdnc_max=', cdnc_max_m3/1.E6
     195  ENDIF
    175196
    176197  d_rei_dt = (rei_max-rei_min)/81.4
    … …  
    204225        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
    205226        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
    206       END DO
    207     END DO
     227      ENDDO
     228    ENDDO
    208229  ELSE ! of IF (iflag_t_glace.EQ.0)
    209230    DO k = 1, klev
    … …  
    222243        xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
    223244        xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
    224       END DO
    225     END DO
     245      ENDDO
     246    ENDDO
    226247  ENDIF
    227248
    … …  
    232253    DO k = 1, klev
    233254      DO i = 1, klon
    234 
    235255        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
    236256        ! Cloud droplet number concentration (CDNC) is restricted
    237257        ! to be within [20, 1000 cm^3]
    238258
    239         ! --present-day case
    240         cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
    241           1.E-4))/log(10.))*1.E6 !-m-3
    242         cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
    243 
    244259        ! --pre-industrial case
    245260        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
    246261          1.E-4))/log(10.))*1.E6 !-m-3
    247         cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
     262        cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
     263
     264      ENDDO
     265    ENDDO
     266
     267    !--flag_aerosol=7 => MACv2SP climatology
     268    !--in this case there is an enhancement factor
     269    IF (flag_aerosol .EQ. 7) THEN
     270
     271      !--present-day
     272      DO k = 1, klev
     273        DO i = 1, klon
     274          cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
     275        ENDDO
     276      ENDDO
     277
     278    !--standard case
     279    ELSE
     280
     281      DO k = 1, klev
     282        DO i = 1, klon
     283
     284          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
     285          ! Cloud droplet number concentration (CDNC) is restricted
     286          ! to be within [20, 1000 cm^3]
     287
     288          ! --present-day case
     289          cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
     290            1.E-4))/log(10.))*1.E6 !-m-3
     291          cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
     292
     293        ENDDO
     294      ENDDO
     295
     296    ENDIF !--flag_aerosol
     297
     298    !--computing cloud droplet size
     299    DO k = 1, klev
     300      DO i = 1, klon
    248301
    249302        ! --present-day case
    … …  
    280333            zfiwp_var*(3.448E-03+2.431/rei)
    281334
    282         END IF
    283 
    284       END DO
    285     END DO
     335        ENDIF
     336
     337      ENDDO
     338    ENDDO
    286339
    287340  ELSE !--not ok_cdnc
    … …  
    293346        rad_chaud(i, k) = rad_chau2
    294347        rad_chaud_pi(i, k) = rad_chau2
    295       END DO
    296     END DO
     348      ENDDO
     349    ENDDO
    297350    DO k = min(3, klev) + 1, klev
    298351      DO i = 1, klon
    299352        rad_chaud(i, k) = rad_chau1
    300353        rad_chaud_pi(i, k) = rad_chau1
    301       END DO
    302     END DO
    303 
    304   END IF !--ok_cdnc
     354      ENDDO
     355    ENDDO
     356
     357  ENDIF !--ok_cdnc
    305358
    306359  ! --computation of cloud optical depth and emissivity
    … …  
    377430        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
    378431
    379       END IF
     432      ENDIF
    380433
    381434      reice(i, k) = rei
    … …  
    384437      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
    385438
    386     END DO
    387   END DO
     439    ENDDO
     440  ENDDO
    388441
    389442  ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
    … …  
    394447        pcldtaupi(i, k) = pcltau(i, k)
    395448        reice_pi(i, k) = reice(i, k)
    396       END DO
    397     END DO
    398   END IF
     449      ENDDO
     450    ENDDO
     451  ENDIF
    399452
    400453  DO k = 1, klev
    … …  
    403456      reliq_pi(i, k) = rad_chaud_pi(i, k)
    404457      reice_pi(i, k) = reice(i, k)
    405     END DO
    406   END DO
     458    ENDDO
     459  ENDDO
    407460
    408461  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
    … …  
    420473    pcl(i) = 1.0
    421474    pctlwp(i) = 0.0
    422   END DO
     475  ENDDO
    423476
    424477  ! --calculation of liquid water path
    … …  
    427480    DO i = 1, klon
    428481      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k)
    429     END DO
    430   END DO
     482    ENDDO
     483  ENDDO
    431484
    432485  ! --calculation of cloud properties with cloud overlap
    … …  
    450503            (i),kind=8),1.-zepsec))
    451504          zcloudl(i) = pclc(i, k)
    452         END IF
     505        ENDIF
    453506        zcloud(i) = pclc(i, k)
    454       END DO
    455     END DO
     507      ENDDO
     508    ENDDO
    456509  ELSE IF (novlp==2) THEN
    457510    DO k = klev, 1, -1
    … …  
    465518        ELSE IF (paprs(i,k)>=prlmc) THEN
    466519          pcl(i) = min(pclc(i,k), pcl(i))
    467         END IF
    468       END DO
    469     END DO
     520        ENDIF
     521      ENDDO
     522    ENDDO
    470523  ELSE IF (novlp==3) THEN
    471524    DO k = klev, 1, -1
    … …  
    479532        ELSE IF (paprs(i,k)>=prlmc) THEN
    480533          pcl(i) = pcl(i)*(1.0-pclc(i,k))
    481         END IF
    482       END DO
    483     END DO
    484   END IF
     534        ENDIF
     535      ENDDO
     536    ENDDO
     537  ENDIF
    485538
    486539  DO i = 1, klon
    … …  
    488541    pcm(i) = 1. - pcm(i)
    489542    pcl(i) = 1. - pcl(i)
    490   END DO
     543  ENDDO
    491544
    492545  ! ========================================================
    … …  
    509562        ELSE
    510563          lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
    511         END IF
     564        ENDIF
    512565        scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
    513       END DO
    514     END DO
     566      ENDDO
     567    ENDDO
    515568
    516569    DO i = 1, klon
    … …  
    520573      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
    521574      IF (novlp.EQ.2) tcc(i) = 0.
    522     END DO
     575    ENDDO
    523576
    524577    DO i = 1, klon
    … …  
    534587              WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM'
    535588              first = .FALSE.
    536             END IF
     589            ENDIF
    537590            flag_max = -1.
    538591            ftmp(i) = max(tcc(i), pclc(i,k))
    539           END IF
     592          ENDIF
    540593
    541594          IF (novlp.EQ.3) THEN
    … …  
    543596              WRITE (*, *) 'Hypothese de recouvrement: RANDOM'
    544597              first = .FALSE.
    545             END IF
     598            ENDIF
    546599            flag_max = 1.
    547600            ftmp(i) = tcc(i)*(1-pclc(i,k))
    548           END IF
     601          ENDIF
    549602
    550603          IF (novlp.EQ.1) THEN
    … …  
    554607                &                                          RANDOM'
    555608              first = .FALSE.
    556             END IF
     609            ENDIF
    557610            flag_max = 1.
    558611            ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
    559612              k+1),1.-thres_neb))
    560           END IF
     613          ENDIF
    561614          ! Effective radius of cloud droplet at top of cloud (m)
    562615          reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, &
    … …  
    570623          tcc(i) = ftmp(i)
    571624
    572         END IF ! is there a visible, not-too-small cloud?
    573       END DO ! loop over k
     625        ENDIF ! is there a visible, not-too-small cloud?
     626      ENDDO ! loop over k
    574627
    575628      IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
    576629
    577     END DO ! loop over i
     630    ENDDO ! loop over i
    578631
    579632    ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
    … …  
    586639        lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
    587640        lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
     641        !FC pour la glace (CAUSES)
     642        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*clwcon(i, k) !  glace convective
     643        icc3dstra(i, k)= pclc(i, k)*pqlwp(i, k)*(1-phase3d(i, k))
     644        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
     645        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
     646        !FC (CAUSES)
     647
    588648        ! Compute cloud droplet radius as above in meter
    589649        radius = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3*rpi*1000.* &
    … …  
    596656        reffclws(i, k) = radius
    597657        reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
    598       END DO !klev
    599     END DO !klon
     658      ENDDO !klev
     659    ENDDO !klon
    600660
    601661    ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
    … …  
    609669        lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
    610670        height(i) = height(i) + dh(i, k)
    611       END DO ! klev
     671      ENDDO ! klev
    612672      lcc_integrat(i) = lcc_integrat(i)/height(i)
    613673      IF (lcc_integrat(i)<=1.0E-03) THEN
    … …  
    615675      ELSE
    616676        cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
    617       END IF
    618     END DO ! klon
     677      ENDIF
     678    ENDDO ! klon
    619679
    620680    DO i = 1, klon
    … …  
    626686        IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
    627687        IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
    628       END DO
     688!FC (CAUSES)
     689        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
     690        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
     691!FC (CAUSES)
     692      ENDDO
    629693      IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
    630694      IF (cldncl(i)<=0.0) cldncl(i) = 0.0
    631695      IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
    632696      IF (lcc(i)<=0.0) lcc(i) = 0.0
    633     END DO
    634 
    635   END IF !ok_cdnc
     697    ENDDO
     698
     699  ENDIF !ok_cdnc
     700
     701  first=.false. !to be sure
    636702
    637703  RETURN
Note: See TracChangeset for help on using the changeset viewer.