Sources/phylmd/newmicro.f

module newmicro_m

  IMPLICIT none

contains

  SUBROUTINE newmicro (paprs, play, t, qlwp, clc, cltau, clemi, cldh, &
       cldl, cldm, cldt, ctlwp, flwp, fiwp, flwc, fiwc, ok_aie, sulfate, &
       sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)

    ! From LMDZ4/libf/phylmd/newmicro.F, version 1.2 2004/06/03 09:22:43

    ! Authors: Z. X. Li (LMD/CNRS), Johannes Quaas
    ! Date: 1993/09/10
    ! Objet: calcul de l'épaisseur optique et de l'émissivité des nuages.

    USE conf_phys_m, ONLY: rad_chau1, rad_chau2
    USE dimphy, ONLY: klev, klon
    USE suphec_m, ONLY: rd, rg
    use nr_util, only: pi

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

    REAL, intent(in):: qlwp(:, :) ! (klon, klev)
    ! eau liquide nuageuse dans l'atmosphère (kg/kg)

    REAL, intent(inout):: clc(:, :) ! (klon, klev)
    ! couverture nuageuse pour le rayonnement (0 à 1)

    REAL, intent(out):: cltau(:, :) ! (klon, klev)  épaisseur optique des nuages
    REAL, intent(out):: clemi(:, :) ! (klon, klev) émissivité des nuages (0 à 1)

    REAL, intent(out):: cldh(:), cldl(:), cldm(:), cldt(:) ! (klon)
    REAL, intent(out):: ctlwp(:) ! (klon)
    REAL, intent(out):: flwp(:), fiwp(:) ! (klon)
    REAL, intent(out):: flwc(:, :), fiwc(:, :) ! (klon, klev)
    LOGICAL, intent(in):: ok_aie ! apply aerosol indirect effect

    REAL, intent(in):: sulfate(:, :) ! (klon, klev)
    ! sulfate aerosol mass concentration (micro g m-3)

    REAL, intent(in):: sulfate_pi(:, :) ! (klon, klev)
    ! sulfate aerosol mass concentration (micro g m-3), pre-industrial value

    REAL, intent(in):: bl95_b0, bl95_b1
    ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
    ! B). They link cloud droplet number concentration to aerosol mass
    ! concentration.

    REAL, intent(out):: cldtaupi(:, :) ! (klon, klev)
    ! pre-industrial value of cloud optical thickness, needed for the
    ! diagnosis of the aerosol indirect radiative forcing (see
    ! radlwsw)

    REAL, intent(out):: re(:, :) ! (klon, klev)
    ! cloud droplet effective radius multiplied by fl (micro m)

    REAL, intent(out):: fl(:, :) ! (klon, klev) 
    ! 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:

    REAL, PARAMETER:: cetahb = 0.45, cetamb = 0.8
    INTEGER i, k
    REAL zflwp(klon), fice
    REAL radius, rad_chaud
    REAL, PARAMETER:: coef_chau = 0.13
    REAL, PARAMETER:: seuil_neb = 0.001, t_glace = 273. - 15.
    real rel, tc, rei, zfiwp(klon)
    real k_ice
    real, parameter:: k_ice0 = 0.005 ! units=m2/g
    real, parameter:: DF = 1.66 ! diffusivity factor
    REAL cdnc(klon, klev) ! cloud droplet number concentration (m-3)

    REAL cdnc_pi(klon, klev)
    ! cloud droplet number concentration, pre-industrial value (m-3)

    !-----------------------------------------------------------------

    ! Calculer l'épaisseur optique et l'émissivité des nuages

    loop_horizontal: DO i = 1, klon
       flwp(i) = 0.
       fiwp(i) = 0.

       DO k = 1, klev
          clc(i, k) = MAX(clc(i, k), seuil_neb)

          ! liquid/ice cloud water paths:

          fice = 1. - (t(i, k) - t_glace) / (273.13 - t_glace)
          fice = MIN(MAX(fice, 0.), 1.)

          zflwp(i) = 1000. * (1. - fice) * qlwp(i, k) / clc(i, k) &
               * (paprs(i, k) - paprs(i, k + 1)) / RG
          zfiwp(i) = 1000. * fice * qlwp(i, k) / clc(i, k) &
               * (paprs(i, k) - paprs(i, k + 1)) / RG

          flwp(i) = flwp(i) &
               + (1. - fice) * qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG
          fiwp(i) = fiwp(i) &
               + fice * qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG

          ! Total Liquid/Ice water content
          flwc(i, k) = (1.-fice) * qlwp(i, k)
          fiwc(i, k) = fice * qlwp(i, k)
          ! In-Cloud Liquid/Ice water content

          ! effective cloud droplet radius (microns):

          ! for liquid water clouds: 
          IF (ok_aie) THEN
             cdnc(i, k) = 10.**(bl95_b0 + bl95_b1 &
                  * log10(MAX(sulfate(i, k), 1e-4)) + 6.)
             cdnc_pi(i, k) = 10.**(bl95_b0 + bl95_b1 &
                  * log10(MAX(sulfate_pi(i, k), 1e-4)) + 6.)

             ! Restrict to interval [20, 1000] cm-3:
             cdnc(i, k) = MIN(1000e6, MAX(20e6, cdnc(i, k)))
             cdnc_pi(i, k) = MIN(1000e6, MAX(20e6, cdnc_pi(i, k)))

             ! air density: play(i, k) / (RD * T(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 * ((qlwp(i, k) * play(i, k) / (RD * T(i, k))) &
                  / (4./3. * PI * 1000. * cdnc(i, k)))**(1./3.)

             ! Convert to micro m and set a lower limit:
             rad_chaud = MAX(rad_chaud * 1e6, 5.) 

             ! Pre-industrial cloud optical thickness

             ! "radius" is calculated as rad_chaud above (plus the 
             ! ice cloud contribution) but using cdnc_pi instead of
             ! cdnc.
             radius = 1.1 * ((qlwp(i, k) * play(i, k) / (RD * T(i, k))) &
                  / (4./3. * PI * 1000. * cdnc_pi(i, k)))**(1./3.)
             radius = MAX(radius * 1e6, 5.) 

             tc = t(i, k)-273.15
             rei = merge(3.5, 0.71 * tc + 61.29, tc <= -81.4)
             if (zflwp(i) == 0.) radius = 1. 
             if (zfiwp(i) == 0. .or. rei <= 0.) rei = 1. 
             cldtaupi(i, k) = 3. / 2. * zflwp(i) / radius &
                  + zfiwp(i) * (3.448e-03 + 2.431 / rei)
          else
             rad_chaud = merge(rad_chau2, rad_chau1, k <= 3)
          ENDIF
          ! For output diagnostics

          ! Cloud droplet effective radius (micro m)

          ! 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) = clc(i, k) * (1.-fice) 
          re(i, k) = rad_chaud * fl(i, k)

          rel = rad_chaud
          ! for ice clouds: as a function of the ambiant temperature
          ! (formula used by Iacobellis and Somerville (2000), with an 
          ! asymptotical value of 3.5 microns at T<-81.4 C added to be 
          ! consistent with observations of Heymsfield et al. 1986):
          tc = t(i, k)-273.15
          rei = merge(3.5, 0.71 * tc + 61.29, tc <= -81.4)

          ! cloud optical thickness:

          ! (for liquid clouds, traditional formula, 
          ! for ice clouds, Ebert & Curry (1992)) 

          if (zflwp(i) == 0.) rel = 1. 
          if (zfiwp(i) == 0. .or. rei <= 0.) rei = 1. 
          cltau(i, k) = 3./2. * (zflwp(i)/rel) &
               + zfiwp(i) * (3.448e-03 + 2.431/rei)

          ! cloud infrared emissivity:

          ! (the broadband infrared absorption coefficient is parameterized
          ! as a function of the effective cld droplet radius)

          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
          k_ice = k_ice0 + 1. / rei

          clemi(i, k) = 1. - EXP(- coef_chau * zflwp(i) - DF * k_ice * zfiwp(i))

          if (clc(i, k) <= seuil_neb) then
             clc(i, k) = 0.
             cltau(i, k) = 0.
             clemi(i, k) = 0.
             cldtaupi(i, k) = 0.
          end if

          IF (.NOT. ok_aie) cldtaupi(i, k) = cltau(i, k) 
       ENDDO
    ENDDO loop_horizontal

    ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS

    DO i = 1, klon
       cldt(i)=1.
       cldh(i)=1.
       cldm(i) = 1.
       cldl(i) = 1.
       ctlwp(i) = 0.
    ENDDO

    DO k = klev, 1, -1
       DO i = 1, klon
          ctlwp(i) = ctlwp(i) &
               + qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG
          cldt(i) = cldt(i) * (1.-clc(i, k))
          if (play(i, k) <= cetahb * paprs(i, 1)) &
               cldh(i) = cldh(i) * (1. - clc(i, k))
          if (play(i, k) > cetahb * paprs(i, 1) .AND. &
               play(i, k) <= cetamb * paprs(i, 1)) &
               cldm(i) = cldm(i) * (1.-clc(i, k))
          if (play(i, k) > cetamb * paprs(i, 1)) &
               cldl(i) = cldl(i) * (1. - clc(i, k))
       ENDDO
    ENDDO

    DO i = 1, klon
       cldt(i)=1.-cldt(i)
       cldh(i)=1.-cldh(i)
       cldm(i)=1.-cldm(i)
       cldl(i)=1.-cldl(i)
    ENDDO

  END SUBROUTINE newmicro

end module newmicro_m
1	guez	68	module newmicro_m
2	guez	3
3	guez	52	IMPLICIT none
4	guez	3
5	guez	68	contains
6	guez	3
7	guez	69	SUBROUTINE newmicro (paprs, play, t, qlwp, clc, cltau, clemi, cldh, &
8			cldl, cldm, cldt, ctlwp, flwp, fiwp, flwc, fiwc, ok_aie, sulfate, &
9			sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
10	guez	3
11	guez	68	! From LMDZ4/libf/phylmd/newmicro.F, version 1.2 2004/06/03 09:22:43
12	guez	3
13	guez	69	! Authors: Z. X. Li (LMD/CNRS), Johannes Quaas
14			! Date: 1993/09/10
15			! Objet: calcul de l'épaisseur optique et de l'émissivité des nuages.
16
17			USE conf_phys_m, ONLY: rad_chau1, rad_chau2
18			USE dimphy, ONLY: klev, klon
19			USE suphec_m, ONLY: rd, rg
20			use nr_util, only: pi
21
22			REAL, intent(in):: paprs(:, :) ! (klon, klev+1)
23			real, intent(in):: play(:, :) ! (klon, klev)
24			REAL, intent(in):: t(:, :) ! (klon, klev) temperature
25
26			REAL, intent(in):: qlwp(:, :) ! (klon, klev)
27			! eau liquide nuageuse dans l'atmosphère (kg/kg)
28
29			REAL, intent(inout):: clc(:, :) ! (klon, klev)
30			! couverture nuageuse pour le rayonnement (0 à 1)
31
32			REAL, intent(out):: cltau(:, :) ! (klon, klev) épaisseur optique des nuages
33			REAL, intent(out):: clemi(:, :) ! (klon, klev) émissivité des nuages (0 à 1)
34
35			REAL, intent(out):: cldh(:), cldl(:), cldm(:), cldt(:) ! (klon)
36			REAL, intent(out):: ctlwp(:) ! (klon)
37			REAL, intent(out):: flwp(:), fiwp(:) ! (klon)
38			REAL, intent(out):: flwc(:, :), fiwc(:, :) ! (klon, klev)
39			LOGICAL, intent(in):: ok_aie ! apply aerosol indirect effect
40
41			REAL, intent(in):: sulfate(:, :) ! (klon, klev)
42			! sulfate aerosol mass concentration (micro g m-3)
43
44			REAL, intent(in):: sulfate_pi(:, :) ! (klon, klev)
45			! sulfate aerosol mass concentration (micro g m-3), pre-industrial value
46
47			REAL, intent(in):: bl95_b0, bl95_b1
48			! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
49			! B). They link cloud droplet number concentration to aerosol mass
50			! concentration.
51
52			REAL, intent(out):: cldtaupi(:, :) ! (klon, klev)
53			! pre-industrial value of cloud optical thickness, needed for the
54			! diagnosis of the aerosol indirect radiative forcing (see
55			! radlwsw)
56
57			REAL, intent(out):: re(:, :) ! (klon, klev)
58			! cloud droplet effective radius multiplied by fl (micro m)
59
60			REAL, intent(out):: fl(:, :) ! (klon, klev)
61			! Denominator to re, introduced to avoid problems in the averaging
62			! of the output. fl is the fraction of liquid water clouds within
63			! a grid cell.
64
65			! Local:
66
67			REAL, PARAMETER:: cetahb = 0.45, cetamb = 0.8
68	guez	68	INTEGER i, k
69	guez	69	REAL zflwp(klon), fice
70	guez	68	REAL radius, rad_chaud
71	guez	69	REAL, PARAMETER:: coef_chau = 0.13
72			REAL, PARAMETER:: seuil_neb = 0.001, t_glace = 273. - 15.
73	guez	68	real rel, tc, rei, zfiwp(klon)
74	guez	69	real k_ice
75			real, parameter:: k_ice0 = 0.005 ! units=m2/g
76			real, parameter:: DF = 1.66 ! diffusivity factor
77			REAL cdnc(klon, klev) ! cloud droplet number concentration (m-3)
78	guez	3
79	guez	69	REAL cdnc_pi(klon, klev)
80			! cloud droplet number concentration, pre-industrial value (m-3)
81	guez	3
82	guez	69	!-----------------------------------------------------------------
83	guez	3
84	guez	69	! Calculer l'épaisseur optique et l'émissivité des nuages
85	guez	3
86	guez	69	loop_horizontal: DO i = 1, klon
87			flwp(i) = 0.
88			fiwp(i) = 0.
89
90	guez	68	DO k = 1, klev
91	guez	69	clc(i, k) = MAX(clc(i, k), seuil_neb)
92	guez	3
93	guez	69	! liquid/ice cloud water paths:
94	guez	3
95	guez	69	fice = 1. - (t(i, k) - t_glace) / (273.13 - t_glace)
96			fice = MIN(MAX(fice, 0.), 1.)
97	guez	3
98	guez	69	zflwp(i) = 1000. * (1. - fice) * qlwp(i, k) / clc(i, k) &
99			* (paprs(i, k) - paprs(i, k + 1)) / RG
100			zfiwp(i) = 1000. * fice * qlwp(i, k) / clc(i, k) &
101			* (paprs(i, k) - paprs(i, k + 1)) / RG
102	guez	3
103	guez	69	flwp(i) = flwp(i) &
104			+ (1. - fice) * qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG
105			fiwp(i) = fiwp(i) &
106			+ fice * qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG
107	guez	3
108	guez	69	! Total Liquid/Ice water content
109			flwc(i, k) = (1.-fice) * qlwp(i, k)
110			fiwc(i, k) = fice * qlwp(i, k)
111			! In-Cloud Liquid/Ice water content
112	guez	3
113	guez	69	! effective cloud droplet radius (microns):
114	guez	3
115	guez	69	! for liquid water clouds:
116			IF (ok_aie) THEN
117			cdnc(i, k) = 10.**(bl95_b0 + bl95_b1 &
118	guez	175	* log10(MAX(sulfate(i, k), 1e-4)) + 6.)
119	guez	69	cdnc_pi(i, k) = 10.**(bl95_b0 + bl95_b1 &
120	guez	175	* log10(MAX(sulfate_pi(i, k), 1e-4)) + 6.)
121	guez	52
122	guez	175	! Restrict to interval [20, 1000] cm-3:
123	guez	69	cdnc(i, k) = MIN(1000e6, MAX(20e6, cdnc(i, k)))
124			cdnc_pi(i, k) = MIN(1000e6, MAX(20e6, cdnc_pi(i, k)))
125	guez	52
126	guez	69	! air density: play(i, k) / (RD * T(i, k))
127			! factor 1.1: derive effective radius from volume-mean radius
128			! factor 1000 is the water density
129			! "_chaud" means that this is the CDR for liquid water clouds
130	guez	52
131	guez	69	rad_chaud = 1.1 * ((qlwp(i, k) * play(i, k) / (RD * T(i, k))) &
132			/ (4./3. * PI * 1000. * cdnc(i, k)))**(1./3.)
133	guez	52
134	guez	69	! Convert to micro m and set a lower limit:
135			rad_chaud = MAX(rad_chaud * 1e6, 5.)
136	guez	52
137	guez	69	! Pre-industrial cloud optical thickness
138	guez	52
139	guez	69	! "radius" is calculated as rad_chaud above (plus the
140			! ice cloud contribution) but using cdnc_pi instead of
141			! cdnc.
142			radius = 1.1 * ((qlwp(i, k) * play(i, k) / (RD * T(i, k))) &
143			/ (4./3. * PI * 1000. * cdnc_pi(i, k)))**(1./3.)
144			radius = MAX(radius * 1e6, 5.)
145	guez	52
146	guez	69	tc = t(i, k)-273.15
147			rei = merge(3.5, 0.71 * tc + 61.29, tc <= -81.4)
148			if (zflwp(i) == 0.) radius = 1.
149			if (zfiwp(i) == 0. .or. rei <= 0.) rei = 1.
150			cldtaupi(i, k) = 3. / 2. * zflwp(i) / radius &
151			+ zfiwp(i) * (3.448e-03 + 2.431 / rei)
152			else
153			rad_chaud = merge(rad_chau2, rad_chau1, k <= 3)
154			ENDIF
155			! For output diagnostics
156	guez	52
157	guez	69	! Cloud droplet effective radius (micro m)
158	guez	52
159	guez	69	! we multiply here with f * xl (fraction of liquid water
160			! clouds in the grid cell) to avoid problems in the
161			! averaging of the output.
162			! In the output of IOIPSL, derive the real cloud droplet
163			! effective radius as re/fl
164	guez	52
165	guez	69	fl(i, k) = clc(i, k) * (1.-fice)
166			re(i, k) = rad_chaud * fl(i, k)
167	guez	52
168	guez	69	rel = rad_chaud
169			! for ice clouds: as a function of the ambiant temperature
170			! (formula used by Iacobellis and Somerville (2000), with an
171			! asymptotical value of 3.5 microns at T<-81.4 C added to be
172			! consistent with observations of Heymsfield et al. 1986):
173			tc = t(i, k)-273.15
174			rei = merge(3.5, 0.71 * tc + 61.29, tc <= -81.4)
175	guez	52
176	guez	69	! cloud optical thickness:
177	guez	68
178	guez	69	! (for liquid clouds, traditional formula,
179			! for ice clouds, Ebert & Curry (1992))
180	guez	68
181	guez	69	if (zflwp(i) == 0.) rel = 1.
182			if (zfiwp(i) == 0. .or. rei <= 0.) rei = 1.
183			cltau(i, k) = 3./2. * (zflwp(i)/rel) &
184			+ zfiwp(i) * (3.448e-03 + 2.431/rei)
185	guez	68
186	guez	69	! cloud infrared emissivity:
187
188			! (the broadband infrared absorption coefficient is parameterized
189			! as a function of the effective cld droplet radius)
190
191			! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
192			k_ice = k_ice0 + 1. / rei
193
194			clemi(i, k) = 1. - EXP(- coef_chau * zflwp(i) - DF * k_ice * zfiwp(i))
195
196			if (clc(i, k) <= seuil_neb) then
197			clc(i, k) = 0.
198			cltau(i, k) = 0.
199			clemi(i, k) = 0.
200			cldtaupi(i, k) = 0.
201			end if
202
203			IF (.NOT. ok_aie) cldtaupi(i, k) = cltau(i, k)
204	guez	68	ENDDO
205	guez	69	ENDDO loop_horizontal
206
207	guez	68	! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
208	guez	69
209	guez	68	DO i = 1, klon
210	guez	69	cldt(i)=1.
211			cldh(i)=1.
212			cldm(i) = 1.
213			cldl(i) = 1.
214			ctlwp(i) = 0.
215	guez	68	ENDDO
216	guez	69
217	guez	68	DO k = klev, 1, -1
218			DO i = 1, klon
219	guez	69	ctlwp(i) = ctlwp(i) &
220			+ qlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / RG
221			cldt(i) = cldt(i) * (1.-clc(i, k))
222			if (play(i, k) <= cetahb * paprs(i, 1)) &
223			cldh(i) = cldh(i) * (1. - clc(i, k))
224			if (play(i, k) > cetahb * paprs(i, 1) .AND. &
225			play(i, k) <= cetamb * paprs(i, 1)) &
226			cldm(i) = cldm(i) * (1.-clc(i, k))
227			if (play(i, k) > cetamb * paprs(i, 1)) &
228			cldl(i) = cldl(i) * (1. - clc(i, k))
229	guez	68	ENDDO
230			ENDDO
231	guez	69
232	guez	68	DO i = 1, klon
233	guez	69	cldt(i)=1.-cldt(i)
234			cldh(i)=1.-cldh(i)
235			cldm(i)=1.-cldm(i)
236			cldl(i)=1.-cldl(i)
237	guez	68	ENDDO
238
239			END SUBROUTINE newmicro
240
241			end module newmicro_m