trunk/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	module nuage_m
2
3	IMPLICIT none
4
5	contains
6
7	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
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	use nr_util, only: pi
15	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	! Arguments:
22
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	! ok_aie--input-L-apply aerosol indirect effect or not
44
45	REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3]
46	! sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
47
48	REAL, intent(in):: sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
49	! sulfate_pi-input-R-dito, pre-industrial value
50
51	REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
52	! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
53	! 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	! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
57	! 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	! re------output-R-Cloud droplet effective radius multiplied by fl [um]
62
63	REAL fl(klon, klev) ! xliq * rneb (denominator to re ; fraction of liquid water clouds within the grid cell)
64	! fl------output-R-Denominator to re, introduced to avoid problems in
65	! the averaging of the output. fl is the fraction of liquid
66	! water clouds within a grid cell
67
68	! Local:
69
70	LOGICAL lo
71
72	REAL cetahb, cetamb
73	PARAMETER (cetahb = 0.45, cetamb = 0.80)
74
75	INTEGER i, k
76	REAL zflwp, zfice
77
78	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	!cc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
81	! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
82	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
89	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
92	!--------------------------------------------------------------------
93
94	! Calculer l'epaisseur optique et l'emmissivite des nuages
95
96	DO k = 1, klev
97	DO i = 1, klon
98	rad_chaud = rad_chau1
99	IF (k <= 3) rad_chaud = rad_chau2
100
101	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	zfice = zfice**nexpo
107
108	IF (ok_aie) THEN
109	! Formula "D" of Boucher and Lohmann, Tellus, 1995
110
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	! Cloud droplet number concentration (CDNC) is restricted
114	! to be within [20, 1000 cm^3]
115
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	! 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
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	! Convert to um. CDR shall be at least 3 um.
131
132	rad_chaud = MAX(rad_chaud*1.e6, 3.)
133
134	! For output diagnostics
135
136	! Cloud droplet effective radius [um]
137
138	! 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
144	fl(i, k) = pclc(i, k)*(1.-zfice)
145	re(i, k) = rad_chaud*fl(i, k)
146
147	! Pre-industrial cloud opt thickness
148
149	! "radius" is calculated as rad_chaud above (plus the
150	! ice cloud contribution) but using cdnc_pi instead of
151	! cdnc.
152	radius = MAX(1.1e6 * ((pqlwp(i, k)pplay(i, k)/(RDT(i, k))) &
153	/ (4./3.PI1000.cdnc_pi(i, k)))*(1./3.), &
154	3.) * (1.-zfice) + rad_froid * zfice
155	cldtaupi(i, k) = 3.0/2.0 * zflwp / radius
156
157	END IF ! ok_aie
158
159	radius = rad_chaud * (1.-zfice) + rad_froid * zfice
160	coef = coef_chau * (1.-zfice) + coef_froi * zfice
161	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
168	IF (.NOT.ok_aie) cldtaupi(i, k)=pcltau(i, k)
169	END DO
170	END DO
171
172	! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
173
174	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
182	DO k = klev, 1, -1
183	DO i = 1, klon
184	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	END DO
195	END DO
196
197	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
204	END SUBROUTINE nuage
205
206	end module nuage_m