Sources/phylmd/nuage.f

module nuage_m

  IMPLICIT none

contains

  SUBROUTINE nuage (paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, &
       pcm, pct, pctlwp, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, &
       cldtaupi, re, fl)

    ! From LMDZ4/libf/phylmd/nuage.F, version 1.1.1.1, 2004/05/19 12:53:07

    use dimphy, only: klon, klev
    use nr_util, only: pi
    use SUPHEC_M, only: rd, rg

    ! Author: Z. X. Li (LMD/CNRS)
    ! Date: 1993/09/10
    ! Objet: Calculer \'epaisseur optique et \'emissivit\'e des nuages

    ! Arguments:

    REAL, intent(in):: paprs(klon, klev+1)
    real, intent(in):: pplay(klon, klev)
    REAL, intent(in):: t(klon, klev) ! temperature

    REAL, intent(in):: pqlwp(klon, klev)
    ! eau liquide nuageuse dans l'atmosphere (kg/kg)

    REAL, intent(inout):: pclc(klon, klev)
    ! couverture nuageuse pour le rayonnement (0 \`a 1)

    REAL pcltau(klon, klev), pclemi(klon, klev)
    ! pcltau--output-R-epaisseur optique des nuages
    ! pclemi--output-R-emissivite des nuages (0 a 1)

    REAL pch(klon), pcl(klon), pcm(klon), pct(klon), pctlwp(klon)

    !jq for the aerosol indirect effect
    !jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
    !jq
    LOGICAL ok_aie ! Apply AIE or not?
    ! ok_aie--input-L-apply aerosol indirect effect or not

    REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3]
    ! sulfate-input-R-sulfate aerosol mass concentration [um/m^3]

    REAL, intent(in):: sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
    ! sulfate_pi-input-R-dito, pre-industrial value

    REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
    ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
    ! bl95_b1-input-R-a parameter, may be varied for tests (-"-)

    REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
    ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
    ! needed for the diagnostics of the aerosol indirect
    ! radiative forcing (see radlwsw)

    REAL re(klon, klev) ! cloud droplet effective radius [um]
    ! re------output-R-Cloud droplet effective radius multiplied by fl [um]

    REAL fl(klon, klev) ! xliq * rneb (denominator to re ; fraction of liquid water clouds within the grid cell)
    ! fl------output-R-Denominator to re, introduced to avoid problems in
    ! the averaging of the output. fl is the fraction of liquid
    ! water clouds within a grid cell

    ! Local:

    LOGICAL lo

    REAL cetahb, cetamb
    PARAMETER (cetahb = 0.45, cetamb = 0.80)

    INTEGER i, k
    REAL zflwp, zfice

    REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
    PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
    !cc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
    ! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
    REAL coef, coef_froi, coef_chau
    PARAMETER (coef_chau=0.13, coef_froi=0.09)
    REAL seuil_neb, t_glace
    PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
    INTEGER nexpo ! exponentiel pour glace/eau
    PARAMETER (nexpo=6)

    REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
    REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)

    !--------------------------------------------------------------------
    
    ! Calculer l'epaisseur optique et l'emmissivite des nuages

    DO k = 1, klev
       DO i = 1, klon
          rad_chaud = rad_chau1
          IF (k <= 3) rad_chaud = rad_chau2

          pclc(i, k) = MAX(pclc(i, k), seuil_neb)
          zflwp = 1000.*pqlwp(i, k)/RG/pclc(i, k) &
               *(paprs(i, k)-paprs(i, k+1))
          zfice = 1.0 - (t(i, k)-t_glace) / (273.13-t_glace)
          zfice = MIN(MAX(zfice, 0.0), 1.0)
          zfice = zfice**nexpo

          IF (ok_aie) THEN
             ! Formula "D" of Boucher and Lohmann, Tellus, 1995

             cdnc(i, k) = 10.**(bl95_b0+bl95_b1* &
                  log(MAX(sulfate(i, k), 1.e-4))/log(10.))*1.e6 !-m-3
             ! Cloud droplet number concentration (CDNC) is restricted
             ! to be within [20, 1000 cm^3]

             cdnc(i, k)=MIN(1000.e6, MAX(20.e6, cdnc(i, k)))
             cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1* &
                  log(MAX(sulfate_pi(i, k), 1.e-4))/log(10.))*1.e6 !-m-3
             cdnc_pi(i, k)=MIN(1000.e6, MAX(20.e6, cdnc_pi(i, k)))

             ! air density: pplay(i, k) / (RD * zT(i, k))
             ! factor 1.1: derive effective radius from volume-mean radius
             ! factor 1000 is the water density
             ! _chaud means that this is the CDR for liquid water clouds

             rad_chaud = &
                  1.1 * ((pqlwp(i, k) * pplay(i, k) / (RD * T(i, k))) &
                  / (4./3. * PI * 1000. * cdnc(i, k)))**(1./3.)

             ! Convert to um. CDR shall be at least 3 um.

             rad_chaud = MAX(rad_chaud*1.e6, 3.)

             ! For output diagnostics

             ! Cloud droplet effective radius [um]

             ! we multiply here with f * xl (fraction of liquid water
             ! clouds in the grid cell) to avoid problems in the
             ! averaging of the output.
             ! In the output of IOIPSL, derive the real cloud droplet
             ! effective radius as re/fl

             fl(i, k) = pclc(i, k)*(1.-zfice)
             re(i, k) = rad_chaud*fl(i, k)

             ! Pre-industrial cloud opt thickness

             ! "radius" is calculated as rad_chaud above (plus the
             ! ice cloud contribution) but using cdnc_pi instead of
             ! cdnc.
             radius = MAX(1.1e6 * ((pqlwp(i, k)*pplay(i, k)/(RD*T(i, k))) &
                  / (4./3.*PI*1000.*cdnc_pi(i, k)))**(1./3.), &
                  3.) * (1.-zfice) + rad_froid * zfice
             cldtaupi(i, k) = 3.0/2.0 * zflwp / radius

          END IF ! ok_aie

          radius = rad_chaud * (1.-zfice) + rad_froid * zfice
          coef = coef_chau * (1.-zfice) + coef_froi * zfice
          pcltau(i, k) = 3.0/2.0 * zflwp / radius
          pclemi(i, k) = 1.0 - EXP(- coef * zflwp)
          lo = (pclc(i, k) <= seuil_neb)
          IF (lo) pclc(i, k) = 0.0
          IF (lo) pcltau(i, k) = 0.0
          IF (lo) pclemi(i, k) = 0.0

          IF (.NOT.ok_aie) cldtaupi(i, k)=pcltau(i, k)
       END DO
    END DO

    ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS

    DO i = 1, klon
       pct(i)=1.0
       pch(i)=1.0
       pcm(i) = 1.0
       pcl(i) = 1.0
       pctlwp(i) = 0.0
    END DO

    DO k = klev, 1, -1
       DO i = 1, klon
          pctlwp(i) = pctlwp(i) &
               + pqlwp(i, k)*(paprs(i, k)-paprs(i, k+1))/RG
          pct(i) = pct(i)*(1.0-pclc(i, k))
          if (pplay(i, k) <= cetahb*paprs(i, 1)) &
               pch(i) = pch(i)*(1.0-pclc(i, k))
          if (pplay(i, k) > cetahb*paprs(i, 1) .AND. &
               pplay(i, k) <= cetamb*paprs(i, 1)) &
               pcm(i) = pcm(i)*(1.0-pclc(i, k))
          if (pplay(i, k) > cetamb*paprs(i, 1)) &
               pcl(i) = pcl(i)*(1.0-pclc(i, k))
       END DO
    END DO

    DO i = 1, klon
       pct(i)=1.-pct(i)
       pch(i)=1.-pch(i)
       pcm(i)=1.-pcm(i)
       pcl(i)=1.-pcl(i)
    END DO

  END SUBROUTINE nuage

end module nuage_m
1	guez	175	module nuage_m
2
3	guez	51	IMPLICIT none
4	guez	3
5	guez	175	contains
6	guez	3
7	guez	175	SUBROUTINE nuage (paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, &
8			pcm, pct, pctlwp, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, &
9			cldtaupi, re, fl)
10	guez	213
11			! From LMDZ4/libf/phylmd/nuage.F, version 1.1.1.1, 2004/05/19 12:53:07
12
13			use dimphy, only: klon, klev
14	guez	175	use nr_util, only: pi
15	guez	213	use SUPHEC_M, only: rd, rg
16
17			! Author: Z. X. Li (LMD/CNRS)
18			! Date: 1993/09/10
19			! Objet: Calculer \'epaisseur optique et \'emissivit\'e des nuages
20
21	guez	175	! Arguments:
22	guez	213
23			REAL, intent(in):: paprs(klon, klev+1)
24			real, intent(in):: pplay(klon, klev)
25			REAL, intent(in):: t(klon, klev) ! temperature
26
27			REAL, intent(in):: pqlwp(klon, klev)
28			! eau liquide nuageuse dans l'atmosphere (kg/kg)
29
30			REAL, intent(inout):: pclc(klon, klev)
31			! couverture nuageuse pour le rayonnement (0 \`a 1)
32
33			REAL pcltau(klon, klev), pclemi(klon, klev)
34			! pcltau--output-R-epaisseur optique des nuages
35			! pclemi--output-R-emissivite des nuages (0 a 1)
36
37			REAL pch(klon), pcl(klon), pcm(klon), pct(klon), pctlwp(klon)
38
39			!jq for the aerosol indirect effect
40			!jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
41			!jq
42			LOGICAL ok_aie ! Apply AIE or not?
43	guez	175	! ok_aie--input-L-apply aerosol indirect effect or not
44	guez	213
45			REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3]
46	guez	175	! sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
47	guez	213
48			REAL, intent(in):: sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
49	guez	175	! sulfate_pi-input-R-dito, pre-industrial value
50	guez	213
51			REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
52	guez	175	! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
53	guez	213	! bl95_b1-input-R-a parameter, may be varied for tests (-"-)
54
55			REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
56	guez	175	! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
57	guez	213	! needed for the diagnostics of the aerosol indirect
58			! radiative forcing (see radlwsw)
59
60			REAL re(klon, klev) ! cloud droplet effective radius [um]
61	guez	175	! re------output-R-Cloud droplet effective radius multiplied by fl [um]
62	guez	213
63			REAL fl(klon, klev) ! xliq * rneb (denominator to re ; fraction of liquid water clouds within the grid cell)
64	guez	175	! fl------output-R-Denominator to re, introduced to avoid problems in
65	guez	213	! the averaging of the output. fl is the fraction of liquid
66			! water clouds within a grid cell
67
68			! Local:
69
70	guez	175	LOGICAL lo
71	guez	213
72	guez	175	REAL cetahb, cetamb
73			PARAMETER (cetahb = 0.45, cetamb = 0.80)
74	guez	213
75	guez	175	INTEGER i, k
76			REAL zflwp, zfice
77	guez	213
78	guez	175	REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
79			PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
80	guez	213	!cc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
81			! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
82	guez	175	REAL coef, coef_froi, coef_chau
83			PARAMETER (coef_chau=0.13, coef_froi=0.09)
84			REAL seuil_neb, t_glace
85			PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
86			INTEGER nexpo ! exponentiel pour glace/eau
87			PARAMETER (nexpo=6)
88	guez	51
89	guez	213	REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
90			REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
91	guez	51
92	guez	213	!--------------------------------------------------------------------
93
94			! Calculer l'epaisseur optique et l'emmissivite des nuages
95	guez	51
96	guez	175	DO k = 1, klev
97			DO i = 1, klon
98			rad_chaud = rad_chau1
99	guez	213	IF (k <= 3) rad_chaud = rad_chau2
100	guez	51
101	guez	213	pclc(i, k) = MAX(pclc(i, k), seuil_neb)
102			zflwp = 1000.*pqlwp(i, k)/RG/pclc(i, k) &
103			*(paprs(i, k)-paprs(i, k+1))
104			zfice = 1.0 - (t(i, k)-t_glace) / (273.13-t_glace)
105			zfice = MIN(MAX(zfice, 0.0), 1.0)
106	guez	175	zfice = zfice**nexpo
107	guez	51
108	guez	175	IF (ok_aie) THEN
109			! Formula "D" of Boucher and Lohmann, Tellus, 1995
110	guez	213
111			cdnc(i, k) = 10.*(bl95_b0+bl95_b1 &
112			log(MAX(sulfate(i, k), 1.e-4))/log(10.))*1.e6 !-m-3
113	guez	175	! Cloud droplet number concentration (CDNC) is restricted
114			! to be within [20, 1000 cm^3]
115	guez	213
116			cdnc(i, k)=MIN(1000.e6, MAX(20.e6, cdnc(i, k)))
117			cdnc_pi(i, k) = 10.*(bl95_b0+bl95_b1 &
118			log(MAX(sulfate_pi(i, k), 1.e-4))/log(10.))*1.e6 !-m-3
119			cdnc_pi(i, k)=MIN(1000.e6, MAX(20.e6, cdnc_pi(i, k)))
120
121			! air density: pplay(i, k) / (RD * zT(i, k))
122	guez	175	! factor 1.1: derive effective radius from volume-mean radius
123			! factor 1000 is the water density
124			! _chaud means that this is the CDR for liquid water clouds
125	guez	213
126			rad_chaud = &
127			1.1 * ((pqlwp(i, k) * pplay(i, k) / (RD * T(i, k))) &
128			/ (4./3. * PI * 1000. * cdnc(i, k)))**(1./3.)
129
130	guez	175	! Convert to um. CDR shall be at least 3 um.
131	guez	213
132	guez	175	rad_chaud = MAX(rad_chaud*1.e6, 3.)
133	guez	51
134	guez	175	! For output diagnostics
135	guez	213
136	guez	175	! Cloud droplet effective radius [um]
137	guez	213
138	guez	175	! we multiply here with f * xl (fraction of liquid water
139			! clouds in the grid cell) to avoid problems in the
140			! averaging of the output.
141			! In the output of IOIPSL, derive the real cloud droplet
142			! effective radius as re/fl
143	guez	51
144	guez	213	fl(i, k) = pclc(i, k)*(1.-zfice)
145			re(i, k) = rad_chaud*fl(i, k)
146
147	guez	175	! Pre-industrial cloud opt thickness
148	guez	213
149	guez	175	! "radius" is calculated as rad_chaud above (plus the
150			! ice cloud contribution) but using cdnc_pi instead of
151			! cdnc.
152	guez	213	radius = MAX(1.1e6 * ((pqlwp(i, k)pplay(i, k)/(RDT(i, k))) &
153			/ (4./3.PI1000.cdnc_pi(i, k)))*(1./3.), &
154	guez	175	3.) * (1.-zfice) + rad_froid * zfice
155	guez	213	cldtaupi(i, k) = 3.0/2.0 * zflwp / radius
156	guez	51
157	guez	213	END IF ! ok_aie
158	guez	175
159			radius = rad_chaud * (1.-zfice) + rad_froid * zfice
160			coef = coef_chau * (1.-zfice) + coef_froi * zfice
161	guez	213	pcltau(i, k) = 3.0/2.0 * zflwp / radius
162			pclemi(i, k) = 1.0 - EXP(- coef * zflwp)
163			lo = (pclc(i, k) <= seuil_neb)
164			IF (lo) pclc(i, k) = 0.0
165			IF (lo) pcltau(i, k) = 0.0
166			IF (lo) pclemi(i, k) = 0.0
167	guez	175
168	guez	213	IF (.NOT.ok_aie) cldtaupi(i, k)=pcltau(i, k)
169	guez	175	END DO
170			END DO
171	guez	213
172	guez	175	! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
173	guez	213
174	guez	175	DO i = 1, klon
175			pct(i)=1.0
176			pch(i)=1.0
177			pcm(i) = 1.0
178			pcl(i) = 1.0
179			pctlwp(i) = 0.0
180			END DO
181	guez	213
182	guez	175	DO k = klev, 1, -1
183			DO i = 1, klon
184	guez	213	pctlwp(i) = pctlwp(i) &
185			+ pqlwp(i, k)*(paprs(i, k)-paprs(i, k+1))/RG
186			pct(i) = pct(i)*(1.0-pclc(i, k))
187			if (pplay(i, k) <= cetahb*paprs(i, 1)) &
188			pch(i) = pch(i)*(1.0-pclc(i, k))
189			if (pplay(i, k) > cetahb*paprs(i, 1) .AND. &
190			pplay(i, k) <= cetamb*paprs(i, 1)) &
191			pcm(i) = pcm(i)*(1.0-pclc(i, k))
192			if (pplay(i, k) > cetamb*paprs(i, 1)) &
193			pcl(i) = pcl(i)*(1.0-pclc(i, k))
194	guez	175	END DO
195			END DO
196	guez	213
197	guez	175	DO i = 1, klon
198			pct(i)=1.-pct(i)
199			pch(i)=1.-pch(i)
200			pcm(i)=1.-pcm(i)
201			pcl(i)=1.-pcl(i)
202			END DO
203	guez	213
204	guez	175	END SUBROUTINE nuage
205
206			end module nuage_m