Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_cloud_optics_prop.f90 @ 5319

Last change on this file since 5319 was 5268, checked in by abarral, 2 weeks ago
.f90 <-> .F90 depending on cpp key use
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.1 KB

Rev	Line
[1279]	1	! $Id: lmdz_cloud_optics_prop.f90 5268 2024-10-23 17:02:39Z fairhead $
[4683]	2	MODULE lmdz_cloud_optics_prop
[1523]	3
[4683]	4	CONTAINS
	5
[4704]	6	SUBROUTINE cloud_optics_prop(klon, klev, paprs, pplay, temp, radocond, picefra, pclc, &
	7	pcltau, pclemi, pch, pcl, pcm, pct, radocondwp, xflwp, xfiwp, xflwc, xfiwc, &
[4692]	8	mass_solu_aero, mass_solu_aero_pi, pcldtaupi, distcltop, temp_cltop, re, fl, reliq, reice, &
	9	reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
	10	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
[4704]	11	icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon)
[524]	12
[4692]	13	USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
	14	USE lmdz_cloud_optics_prop_ini , ONLY : lunout
	15	USE lmdz_cloud_optics_prop_ini , ONLY : bl95_b0, bl95_b1
	16	USE lmdz_cloud_optics_prop_ini , ONLY : latitude_deg
	17	USE lmdz_cloud_optics_prop_ini , ONLY : iflag_t_glace
	18	USE lmdz_cloud_optics_prop_ini , ONLY : cdnc_max, cdnc_max_m3
	19	USE lmdz_cloud_optics_prop_ini , ONLY : cdnc_min, cdnc_min_m3
[4707]	20	USE lmdz_cloud_optics_prop_ini , ONLY : thres_tau, thres_neb
	21	USE lmdz_cloud_optics_prop_ini , ONLY : prmhc, prlmc
	22	USE lmdz_cloud_optics_prop_ini , ONLY : coef_froi, coef_chau
	23	USE lmdz_cloud_optics_prop_ini , ONLY : seuil_neb
	24	USE lmdz_cloud_optics_prop_ini , ONLY : t_glace_min_old, t_glace_max_old
	25	USE lmdz_cloud_optics_prop_ini , ONLY : k_ice0, df
	26	USE lmdz_cloud_optics_prop_ini , ONLY : rg, rd, rpi
	27	USE lmdz_cloud_optics_prop_ini , ONLY : rad_chau1, rad_chau2, rad_froid, iflag_rei
	28	USE lmdz_cloud_optics_prop_ini , ONLY : ok_icefra_lscp, rei_max, rei_min
	29	USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
	30
[3245]	31
[3265]	32
[3999]	33
[1992]	34	IMPLICIT NONE
	35	! ======================================================================
[4704]	36	! Authors: Z.X. Li (LMD/CNRS) date: 19930910
	37	! O.Boucher (LMD/CNRS) mise a jour en 201212
	38	! I. Musat (LMD/CNRS) : prise en compte de la meme hypothese
	39	! de recouvrement pour les nuages que pour
	40	! le rayonnement rrtm via le parametre
	41	! novlp de radopt.h : 20160721
	42	! L.Fairheard, E.Vignon, JB Madeleine, L. Raillard, A. Idelkadi
	43	! M. Coulon-Decorzens: replayisation of the routine + cleaning
	44	! and commentaries
	45	!
	46	! Aim: compute condensate optical properties,
	47	! cloud optical depth and emissivity
[1992]	48	! ======================================================================
[4704]	49
	50	! List of arguments
	51	!------------------
[1992]	52
[4704]	53	! input:
	54	INTEGER, INTENT(IN) :: klon, klev ! number of horizontal and vertical grid points
	55	REAL, INTENT(IN) :: paprs(klon, klev+1)! pressure at bottom interfaces [Pa]
	56	REAL, INTENT(IN) :: pplay(klon, klev) ! pressure at the middle of layers [Pa]
	57	REAL, INTENT(IN) :: temp(klon, klev) ! temperature [K]
	58	REAL, INTENT(IN) :: radocond(klon, klev) ! cloud condensed water seen by radiation [kg/kg]
	59	REAL, INTENT(IN) :: picefra(klon,klev) ! ice fraction in clouds from large scale condensation scheme [-]
	60	REAL, INTENT(IN) :: rnebcon(klon,klev) ! convection cloud fraction [-]
	61	REAL, INTENT(IN) :: ccwcon(klon,klev) ! condensed water from deep convection [kg/kg]
	62	! jq for the aerosol indirect effect
	63	! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
	64	REAL, INTENT(IN) :: mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3]
	65	REAL, INTENT(IN) :: mass_solu_aero_pi(klon, klev) ! - (pre-industrial value)
	66	REAL, INTENT(IN) :: dNovrN(klon) ! enhancement factor for cdnc
	67	REAL, INTENT(OUT) :: distcltop(klon,klev) ! distance from large scale cloud top [m]
	68	REAL, INTENT(OUT) :: temp_cltop(klon,klev)!temperature at large scale cloud top [K]
[1992]	69
[4704]	70	LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
[1992]	71
[4704]	72	! inout:
	73	REAL, INTENT(INOUT) :: pclc(klon, klev) ! cloud fraction for radiation [-]
[1992]	74
[4704]	75	! out:
	76	REAL, INTENT(OUT) :: pct(klon) ! 2D total cloud cover [-]
	77	REAL, INTENT(OUT) :: pcl(klon) ! 2D low cloud cover [-]
	78	REAL, INTENT(OUT) :: pcm(klon) ! 2D mid cloud cover [-]
	79	REAL, INTENT(OUT) :: pch(klon) ! 2D high cloud cover [-]
	80	REAL, INTENT(OUT) :: radocondwp(klon) ! total condensed water path (seen by radiation) [kg/m2]
	81	REAL, INTENT(OUT) :: xflwp(klon) ! liquid water path (seen by radiation) [kg/m2]
	82	REAL, INTENT(OUT) :: xfiwp(klon) ! ice water path (seen by radiation) [kg/m2]
	83	REAL, INTENT(OUT) :: xflwc(klon, klev) ! liquid water content seen by radiation [kg/kg]
	84	REAL, INTENT(OUT) :: xfiwc(klon, klev) ! ice water content seen by radiation [kg/kg]
	85	REAL, INTENT(OUT) :: re(klon, klev) ! cloud droplet effective radius multiplied by fl
	86	REAL, INTENT(OUT) :: fl(klon, klev) ! xliq * rneb, denominator to re; fraction of liquid water clouds
	87	! introduced to avoid problems in the averaging of the output
	88	! water clouds within a grid cell
[4692]	89
[4704]	90	REAL, INTENT(OUT) :: pcltau(klon, klev) ! cloud optical depth [m]
	91	REAL, INTENT(OUT) :: pclemi(klon, klev) ! cloud emissivity [-]
	92	REAL, INTENT(OUT) :: pcldtaupi(klon, klev) ! pre-industrial value of cloud optical thickness, ie.
	93	! values of optical thickness that does not account
	94	! for aerosol effects on cloud droplet radius [m]
[1992]	95
[4704]	96	REAL, INTENT(OUT) :: reliq(klon, klev) ! liquid droplet effective radius [m]
	97	REAL, INTENT(OUT) :: reice(klon, klev) ! ice effective radius [m]
	98	REAL, INTENT(OUT) :: reliq_pi(klon, klev)! liquid droplet effective radius [m], pre-industrial
	99	REAL, INTENT(OUT) :: reice_pi(klon, klev)! ice effective radius [m], pre-industrial
	100	REAL, INTENT(OUT) :: scdnc(klon, klev) ! cloud droplet number concentration, mean over the whole mesh [m-3]
	101	REAL, INTENT(OUT) :: cldncl(klon) ! cloud droplet number concentration at top of cloud [m-3]
	102	REAL, INTENT(OUT) :: reffclwtop(klon) ! effective radius of cloud droplet at top of cloud [m]
	103	REAL, INTENT(OUT) :: lcc(klon) ! liquid Cloud Content at top of cloud [kg/kg]
	104	REAL, INTENT(OUT) :: reffclws(klon, klev)! stratiform cloud droplet effective radius
	105	REAL, INTENT(OUT) :: reffclwc(klon, klev)! convective cloud droplet effective radius
	106	REAL, INTENT(OUT) :: cldnvi(klon) ! column Integrated cloud droplet Number [/m2]
	107	REAL, INTENT(OUT) :: lcc3d(klon, klev) ! cloud fraction for liquid part only [-]
	108	REAL, INTENT(OUT) :: lcc3dcon(klon, klev)! cloud fraction for liquid part only, convective clouds [-]
	109	REAL, INTENT(OUT) :: lcc3dstra(klon, klev)!cloud fraction for liquid part only, stratiform clouds [-]
	110	REAL, INTENT(OUT) :: icc3dcon(klon, klev)! cloud fraction for liquid part only, convective clouds [-]
	111	REAL, INTENT(OUT) :: icc3dstra(klon, klev)! cloud fraction for ice part only, stratiform clouds [-]
[4715]	112	REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
[1992]	113
[4704]	114	! Local variables
	115	!----------------
	116
[1992]	117	LOGICAL, SAVE :: first = .TRUE.
	118	!$OMP THREADPRIVATE(FIRST)
	119	INTEGER flag_max
	120
	121	! threshold PARAMETERs
	122	REAL phase3d(klon, klev)
	123	REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
	124	LOGICAL lo
	125	INTEGER i, k
	126	REAL radius
	127
	128
[4114]	129	REAL rel, tc, rei, iwc, dei, deimin, deimax
[4707]	130	REAL k_ice
[1992]	131
	132	! jq for the aerosol indirect effect
	133	! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
	134	REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
	135	REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
	136	REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
	137
	138	! IM cf. CR:parametres supplementaires
[3999]	139	REAL dzfice(klon,klev)
[1992]	140	REAL zclear(klon)
	141	REAL zcloud(klon)
	142	REAL zcloudh(klon)
	143	REAL zcloudm(klon)
	144	REAL zcloudl(klon)
	145	REAL rhodz(klon, klev) !--rho*dz pour la couche
	146	REAL zrho(klon, klev) !--rho pour la couche
	147	REAL dh(klon, klev) !--dz pour la couche
	148	REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
	149	REAL rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
	150	REAL zflwp_var, zfiwp_var
	151	REAL d_rei_dt
	152
	153
	154	! FH : 2011/05/24
	155	! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
	156	! to be used for a temperature in celcius T(°C) < 0
	157	! rei=rei_min for T(°C) < -81.4
	158	! Calcul de la pente de la relation entre rayon effective des cristaux
[4704]	159	! et la température Pour retrouver les résultats numériques de la version d'origine,
[1992]	160	! on impose 0.71 quand on est proche de 0.71
	161	d_rei_dt = (rei_max-rei_min)/81.4
	162	IF (abs(d_rei_dt-0.71)<1.E-4) d_rei_dt = 0.71
	163
	164	! Calculer l'epaisseur optique et l'emmissivite des nuages
	165	! IM inversion des DO
	166
	167	xflwp = 0.D0
	168	xfiwp = 0.D0
	169	xflwc = 0.D0
	170	xfiwc = 0.D0
	171
	172	reliq = 0.
	173	reice = 0.
	174	reliq_pi = 0.
	175	reice_pi = 0.
	176
[2006]	177	IF (iflag_t_glace.EQ.0) THEN
	178	DO k = 1, klev
	179	DO i = 1, klon
	180	! -layer calculation
	181	rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
[4704]	182	zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
[2006]	183	dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
	184	! -Fraction of ice in cloud using a linear transition
[4704]	185	icefrac_optics(i, k) = 1.0 - (temp(i,k)-t_glace_min_old)/(t_glace_max_old-t_glace_min_old)
	186	icefrac_optics(i, k) = min(max(icefrac_optics(i,k),0.0), 1.0)
[2006]	187	! -IM Total Liquid/Ice water content
[4704]	188	xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
	189	xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
[3274]	190	ENDDO
	191	ENDDO
[2006]	192	ELSE ! of IF (iflag_t_glace.EQ.0)
	193	DO k = 1, klev
[2077]	194
[3999]	195
[4715]	196	!!$ IF (ok_new_lscp) THEN
	197	!!$ CALL icefrac_lscp(klon,temp(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),icefrac_optics(:,k),dzfice(:,k))
	198	!!$ ELSE
	199	!!$ CALL icefrac_lsc(klon,temp(:,k),pplay(:,k)/paprs(:,1),icefrac_optics(:,k))
	200	!!$ ENDIF
[3999]	201
[2006]	202	DO i = 1, klon
[3999]	203
	204	IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
	205	! EV: take the ice fraction directly from the lscp code
	206	! consistent only for non convective grid points
	207	! critical for mixed phase clouds
[4704]	208	icefrac_optics(i,k)=picefra(i,k)
[3999]	209	ENDIF
	210
[2006]	211	! -layer calculation
	212	rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2
[4704]	213	zrho(i, k) = pplay(i, k)/temp(i, k)/rd ! kg/m3
[2006]	214	dh(i, k) = rhodz(i, k)/zrho(i, k) ! m
	215	! -IM Total Liquid/Ice water content
[4704]	216	xflwc(i, k) = (1.-icefrac_optics(i,k))*radocond(i, k)
	217	xfiwc(i, k) = icefrac_optics(i, k)*radocond(i, k)
[3274]	218	ENDDO
	219	ENDDO
[2006]	220	ENDIF
[1992]	221
[4704]	222
	223
	224
	225
	226
[1992]	227	IF (ok_cdnc) THEN
	228
	229	! --we compute cloud properties as a function of the aerosol load
	230
	231	DO k = 1, klev
	232	DO i = 1, klon
	233	! Formula "D" of Boucher and Lohmann, Tellus, 1995
	234	! Cloud droplet number concentration (CDNC) is restricted
	235	! to be within [20, 1000 cm^3]
	236
	237	! --pre-industrial case
	238	cdnc_pi(i, k) = 10.*(bl95_b0+bl95_b1log(max(mass_solu_aero_pi(i,k), &
	239	1.E-4))/log(10.))*1.E6 !-m-3
[3281]	240	cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
[1992]	241
[3274]	242	ENDDO
	243	ENDDO
	244
	245	!--flag_aerosol=7 => MACv2SP climatology
	246	!--in this case there is an enhancement factor
	247	IF (flag_aerosol .EQ. 7) THEN
	248
	249	!--present-day
	250	DO k = 1, klev
	251	DO i = 1, klon
	252	cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
	253	ENDDO
	254	ENDDO
	255
	256	!--standard case
	257	ELSE
	258
	259	DO k = 1, klev
	260	DO i = 1, klon
	261
	262	! Formula "D" of Boucher and Lohmann, Tellus, 1995
	263	! Cloud droplet number concentration (CDNC) is restricted
	264	! to be within [20, 1000 cm^3]
	265
	266	! --present-day case
	267	cdnc(i, k) = 10.*(bl95_b0+bl95_b1log(max(mass_solu_aero(i,k), &
	268	1.E-4))/log(10.))*1.E6 !-m-3
[3281]	269	cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
[3274]	270
	271	ENDDO
	272	ENDDO
	273
	274	ENDIF !--flag_aerosol
	275
	276	!--computing cloud droplet size
	277	DO k = 1, klev
	278	DO i = 1, klon
	279
[1992]	280	! --present-day case
[4704]	281	rad_chaud(i, k) = 1.1((radocond(i,k)pplay(i, &
	282	k)/(rdtemp(i,k)))/(4./3rpi1000.cdnc(i,k)))**(1./3.)
[1992]	283	rad_chaud(i, k) = max(rad_chaud(i,k)*1.E6, 5.)
	284
	285	! --pre-industrial case
[4704]	286	rad_chaud_pi(i, k) = 1.1((radocond(i,k)pplay(i, &
	287	k)/(rdtemp(i,k)))/(4./3.rpi1000.cdnc_pi(i,k)))**(1./3.)
[1992]	288	rad_chaud_pi(i, k) = max(rad_chaud_pi(i,k)*1.E6, 5.)
	289
	290	! --pre-industrial case
	291	! --liquid/ice cloud water paths:
	292	IF (pclc(i,k)<=seuil_neb) THEN
	293
	294	pcldtaupi(i, k) = 0.0
	295
	296	ELSE
	297
[4704]	298	zflwp_var = 1000.(1.-icefrac_optics(i,k))radocond(i, k)/pclc(i, k)* &
[1992]	299	rhodz(i, k)
[4704]	300	zfiwp_var = 1000.icefrac_optics(i, k)radocond(i, k)/pclc(i, k)*rhodz(i, k)
[4114]	301	! Calculation of ice cloud effective radius in micron
	302	IF (iflag_rei .EQ. 1) THEN
	303	! when we account for precipitation in the radiation scheme,
	304	! It is recommended to use the rei formula from Sun and Rikkus 1999 with a revision
	305	! from Sun 2001 (as in the IFS model)
[4704]	306	iwc=icefrac_optics(i, k)radocond(i, k)/pclc(i,k)zrho(i,k)*1000. !in cloud ice water content in g/m3
	307	dei=(1.2351+0.0105(temp(i,k)-273.15))(45.8966(iwc*0.2214) + &
	308	& 0.7957(iwc0.2535)(temp(i,k)-83.15))
[4114]	309	!deimax=155.0
	310	!deimin=20.+40cos(abs(latitude_deg(i))/180.RPI)
	311	!Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
	312	deimax=rei_max*2.0
	313	deimin=2.0rei_min+40cos(abs(latitude_deg(i))/180.*RPI)
	314	dei=min(dei,deimax)
	315	dei=max(dei,deimin)
	316	rei=3.sqrt(3.)/8.dei
	317	ELSE
	318	! Default
	319	! for ice clouds: as a function of the ambiant temperature
	320	! [formula used by Iacobellis and Somerville (2000), with an
	321	! asymptotical value of 3.5 microns at T<-81.4 C added to be
	322	! consistent with observations of Heymsfield et al. 1986]:
	323	! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
	324	! rei_max=61.29
[4704]	325	tc = temp(i, k) - 273.15
[4114]	326	rei = d_rei_dt*tc + rei_max
	327	IF (tc<=-81.4) rei = rei_min
	328	ENDIF
[1992]	329
	330	! -- cloud optical thickness :
	331	! [for liquid clouds, traditional formula,
	332	! for ice clouds, Ebert & Curry (1992)]
	333
	334	IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
	335	pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
	336	zfiwp_var*(3.448E-03+2.431/rei)
	337
[3274]	338	ENDIF
[1992]	339
[3274]	340	ENDDO
	341	ENDDO
[1992]	342
	343	ELSE !--not ok_cdnc
	344
	345	! -prescribed cloud droplet radius
	346
	347	DO k = 1, min(3, klev)
	348	DO i = 1, klon
	349	rad_chaud(i, k) = rad_chau2
	350	rad_chaud_pi(i, k) = rad_chau2
[3274]	351	ENDDO
	352	ENDDO
[1992]	353	DO k = min(3, klev) + 1, klev
	354	DO i = 1, klon
	355	rad_chaud(i, k) = rad_chau1
	356	rad_chaud_pi(i, k) = rad_chau1
[3274]	357	ENDDO
	358	ENDDO
[1992]	359
[3274]	360	ENDIF !--ok_cdnc
[1992]	361
	362	! --computation of cloud optical depth and emissivity
	363	! --in the general case
	364
	365	DO k = 1, klev
	366	DO i = 1, klon
	367
	368	IF (pclc(i,k)<=seuil_neb) THEN
	369
	370	! effective cloud droplet radius (microns) for liquid water clouds:
	371	! For output diagnostics cloud droplet effective radius [um]
	372	! we multiply here with f * xl (fraction of liquid water
	373	! clouds in the grid cell) to avoid problems in the averaging of the
	374	! output.
	375	! In the output of IOIPSL, derive the REAL cloud droplet
	376	! effective radius as re/fl
	377
[4704]	378	fl(i, k) = seuil_neb*(1.-icefrac_optics(i,k))
[1992]	379	re(i, k) = rad_chaud(i, k)*fl(i, k)
	380	rel = 0.
	381	rei = 0.
	382	pclc(i, k) = 0.0
	383	pcltau(i, k) = 0.0
	384	pclemi(i, k) = 0.0
	385
	386	ELSE
	387
	388	! -- liquid/ice cloud water paths:
	389
[4704]	390	zflwp_var = 1000.(1.-icefrac_optics(i,k))radocond(i, k)/pclc(i, k)*rhodz(i, k)
	391	zfiwp_var = 1000.icefrac_optics(i, k)radocond(i, k)/pclc(i, k)*rhodz(i, k)
[1992]	392
	393	! effective cloud droplet radius (microns) for liquid water clouds:
	394	! For output diagnostics cloud droplet effective radius [um]
[4704]	395	! we multiply here with f Effective radius of cloud droplet at top of cloud (m)* xl (fraction of liquid water
[1992]	396	! clouds in the grid cell) to avoid problems in the averaging of the
	397	! output.
	398	! In the output of IOIPSL, derive the REAL cloud droplet
	399	! effective radius as re/fl
	400
[4704]	401	fl(i, k) = pclc(i, k)*(1.-icefrac_optics(i,k))
[1992]	402	re(i, k) = rad_chaud(i, k)*fl(i, k)
	403
	404	rel = rad_chaud(i, k)
	405
[4114]	406	! Calculation of ice cloud effective radius in micron
[1992]	407
	408
[4114]	409	IF (iflag_rei .GT. 0) THEN
	410
	411	! when we account for precipitation in the radiation scheme,
	412	! we use the rei formula from Sun and Rikkus 1999 with a revision
	413	! from Sun 2001 (as in the IFS model)
[4704]	414	iwc=icefrac_optics(i, k)radocond(i, k)/pclc(i,k)zrho(i,k)*1000. !in cloud ice water content in g/m3
	415	dei=(1.2351+0.0105(temp(i,k)-273.15))(45.8966(iwc*0.2214) + &
	416	&0.7957(iwc0.2535)(temp(i,k)-83.15))
[4119]	417	!deimax=155.0
	418	!deimin=20.+40cos(abs(latitude_deg(i))/180.RPI)
	419	!Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
	420	deimax=rei_max*2.0
	421	deimin=2.0rei_min+40cos(abs(latitude_deg(i))/180.*RPI)
[4114]	422	dei=min(dei,deimax)
	423	dei=max(dei,deimin)
	424	rei=3.sqrt(3.)/8.dei
	425
	426	ELSE
	427	! Default
	428	! for ice clouds: as a function of the ambiant temperature
	429	! [formula used by Iacobellis and Somerville (2000), with an
	430	! asymptotical value of 3.5 microns at T<-81.4 C added to be
	431	! consistent with observations of Heymsfield et al. 1986]:
	432	! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
	433	! rei_max=61.29
[4704]	434	tc = temp(i, k) - 273.15
[4114]	435	rei = d_rei_dt*tc + rei_max
	436	IF (tc<=-81.4) rei = rei_min
	437	ENDIF
[1992]	438	! -- cloud optical thickness :
	439	! [for liquid clouds, traditional formula,
	440	! for ice clouds, Ebert & Curry (1992)]
	441
	442	IF (zflwp_var==0.) rel = 1.
	443	IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
	444	pcltau(i, k) = 3.0/2.0(zflwp_var/rel) + zfiwp_var(3.448E-03+2.431/ &
	445	rei)
	446
	447	! -- cloud infrared emissivity:
	448	! [the broadband infrared absorption coefficient is PARAMETERized
	449	! as a function of the effective cld droplet radius]
	450	! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
	451
	452	k_ice = k_ice0 + 1.0/rei
	453
	454	pclemi(i, k) = 1.0 - exp(-coef_chauzflwp_var-dfk_ice*zfiwp_var)
	455
[3274]	456	ENDIF
[1992]	457
	458	reice(i, k) = rei
	459
	460	xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
	461	xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
	462
[3274]	463	ENDDO
	464	ENDDO
[1992]	465
	466	! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
	467
	468	IF (.NOT. ok_cdnc) THEN
	469	DO k = 1, klev
	470	DO i = 1, klon
	471	pcldtaupi(i, k) = pcltau(i, k)
	472	reice_pi(i, k) = reice(i, k)
[3274]	473	ENDDO
	474	ENDDO
	475	ENDIF
[1992]	476
	477	DO k = 1, klev
	478	DO i = 1, klon
	479	reliq(i, k) = rad_chaud(i, k)
	480	reliq_pi(i, k) = rad_chaud_pi(i, k)
	481	reice_pi(i, k) = reice(i, k)
[3274]	482	ENDDO
	483	ENDDO
[1992]	484
	485	! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
	486	! IM cf. CR:test: calcul prenant ou non en compte le recouvrement
	487	! initialisations
	488
	489	DO i = 1, klon
	490	zclear(i) = 1.
	491	zcloud(i) = 0.
	492	zcloudh(i) = 0.
	493	zcloudm(i) = 0.
	494	zcloudl(i) = 0.
	495	pch(i) = 1.0
	496	pcm(i) = 1.0
	497	pcl(i) = 1.0
[4704]	498	radocondwp(i) = 0.0
[3274]	499	ENDDO
[1992]	500
	501	! --calculation of liquid water path
	502
	503	DO k = klev, 1, -1
	504	DO i = 1, klon
[4704]	505	radocondwp(i) = radocondwp(i) + radocond(i, k)*rhodz(i, k)
[3274]	506	ENDDO
	507	ENDDO
[1992]	508
	509	! --calculation of cloud properties with cloud overlap
[4704]	510	! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
	511	! !novlp=1: max-random
	512	! !novlp=2: maximum
	513	! !novlp=3: random
[1992]	514
[4704]	515
[1992]	516	IF (novlp==1) THEN
	517	DO k = klev, 1, -1
	518	DO i = 1, klon
	519	zclear(i) = zclear(i)*(1.-max(pclc(i,k),zcloud(i)))/(1.-min(real( &
	520	zcloud(i),kind=8),1.-zepsec))
	521	pct(i) = 1. - zclear(i)
	522	IF (paprs(i,k)<prmhc) THEN
	523	pch(i) = pch(i)*(1.-max(pclc(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
	524	(i),kind=8),1.-zepsec))
	525	zcloudh(i) = pclc(i, k)
	526	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
	527	pcm(i) = pcm(i)*(1.-max(pclc(i,k),zcloudm(i)))/(1.-min(real(zcloudm &
	528	(i),kind=8),1.-zepsec))
	529	zcloudm(i) = pclc(i, k)
	530	ELSE IF (paprs(i,k)>=prlmc) THEN
	531	pcl(i) = pcl(i)*(1.-max(pclc(i,k),zcloudl(i)))/(1.-min(real(zcloudl &
	532	(i),kind=8),1.-zepsec))
	533	zcloudl(i) = pclc(i, k)
[3274]	534	ENDIF
[1992]	535	zcloud(i) = pclc(i, k)
[3274]	536	ENDDO
	537	ENDDO
[1992]	538	ELSE IF (novlp==2) THEN
	539	DO k = klev, 1, -1
	540	DO i = 1, klon
	541	zcloud(i) = max(pclc(i,k), zcloud(i))
	542	pct(i) = zcloud(i)
	543	IF (paprs(i,k)<prmhc) THEN
	544	pch(i) = min(pclc(i,k), pch(i))
	545	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
	546	pcm(i) = min(pclc(i,k), pcm(i))
	547	ELSE IF (paprs(i,k)>=prlmc) THEN
	548	pcl(i) = min(pclc(i,k), pcl(i))
[3274]	549	ENDIF
	550	ENDDO
	551	ENDDO
[1992]	552	ELSE IF (novlp==3) THEN
	553	DO k = klev, 1, -1
	554	DO i = 1, klon
	555	zclear(i) = zclear(i)*(1.-pclc(i,k))
	556	pct(i) = 1 - zclear(i)
	557	IF (paprs(i,k)<prmhc) THEN
	558	pch(i) = pch(i)*(1.0-pclc(i,k))
	559	ELSE IF (paprs(i,k)>=prmhc .AND. paprs(i,k)<prlmc) THEN
	560	pcm(i) = pcm(i)*(1.0-pclc(i,k))
	561	ELSE IF (paprs(i,k)>=prlmc) THEN
	562	pcl(i) = pcl(i)*(1.0-pclc(i,k))
[3274]	563	ENDIF
	564	ENDDO
	565	ENDDO
	566	ENDIF
[1992]	567
	568	DO i = 1, klon
	569	pch(i) = 1. - pch(i)
	570	pcm(i) = 1. - pcm(i)
	571	pcl(i) = 1. - pcl(i)
[3274]	572	ENDDO
[1992]	573
	574	! ========================================================
	575	! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
	576	! ========================================================
	577	! change by Nicolas Yan (LSCE)
	578	! Cloud Droplet Number Concentration (CDNC) : 3D variable
	579	! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
	580	! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
	581	! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
	582	! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
	583
	584	IF (ok_cdnc) THEN
	585
	586	DO k = 1, klev
	587	DO i = 1, klon
[4704]	588	phase3d(i, k) = 1 - icefrac_optics(i, k)
[1992]	589	IF (pclc(i,k)<=seuil_neb) THEN
	590	lcc3d(i, k) = seuil_neb*phase3d(i, k)
	591	ELSE
	592	lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
[3274]	593	ENDIF
[1992]	594	scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
[3274]	595	ENDDO
	596	ENDDO
[1992]	597
	598	DO i = 1, klon
	599	lcc(i) = 0.
	600	reffclwtop(i) = 0.
	601	cldncl(i) = 0.
[2596]	602	IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
	603	IF (novlp.EQ.2) tcc(i) = 0.
[3274]	604	ENDDO
[1992]	605
	606	DO i = 1, klon
	607	DO k = klev - 1, 1, -1 !From TOA down
	608
	609	! Test, if the cloud optical depth exceeds the necessary
	610	! threshold:
	611
	612	IF (pcltau(i,k)>thres_tau .AND. pclc(i,k)>thres_neb) THEN
	613
[2596]	614	IF (novlp.EQ.2) THEN
[1992]	615	IF (first) THEN
	616	WRITE (, ) 'Hypothese de recouvrement: MAXIMUM'
	617	first = .FALSE.
[3274]	618	ENDIF
[1992]	619	flag_max = -1.
	620	ftmp(i) = max(tcc(i), pclc(i,k))
[3274]	621	ENDIF
[1992]	622
[2596]	623	IF (novlp.EQ.3) THEN
[1992]	624	IF (first) THEN
	625	WRITE (, ) 'Hypothese de recouvrement: RANDOM'
	626	first = .FALSE.
[3274]	627	ENDIF
[1992]	628	flag_max = 1.
	629	ftmp(i) = tcc(i)*(1-pclc(i,k))
[3274]	630	ENDIF
[1992]	631
[2596]	632	IF (novlp.EQ.1) THEN
[1992]	633	IF (first) THEN
	634	WRITE (, ) 'Hypothese de recouvrement: MAXIMUM_ &
	635	& &
	636	& RANDOM'
	637	first = .FALSE.
[3274]	638	ENDIF
[1992]	639	flag_max = 1.
	640	ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
	641	k+1),1.-thres_neb))
[3274]	642	ENDIF
[1992]	643	! Effective radius of cloud droplet at top of cloud (m)
	644	reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)1.0E-06phase3d(i, &
	645	k)(tcc(i)-ftmp(i))flag_max
	646	! CDNC at top of cloud (m-3)
	647	cldncl(i) = cldncl(i) + cdnc(i, k)phase3d(i, k)(tcc(i)-ftmp(i))* &
	648	flag_max
	649	! Liquid Cloud Content at top of cloud
	650	lcc(i) = lcc(i) + phase3d(i, k)(tcc(i)-ftmp(i))flag_max
	651	! Total Cloud Content at top of cloud
	652	tcc(i) = ftmp(i)
	653
[3274]	654	ENDIF ! is there a visible, not-too-small cloud?
	655	ENDDO ! loop over k
[1992]	656
[2596]	657	IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
[1992]	658
[3274]	659	ENDDO ! loop over i
[1992]	660
	661	! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
	662	! REFFCLWS)
	663	DO i = 1, klon
[524]	664	DO k = 1, klev
[1992]	665	! Weight to be used for outputs: eau_liquide*couverture nuageuse
[4704]	666	lcc3dcon(i, k) = rnebcon(i, k)phase3d(i, k)ccwcon(i, k) ! eau liquide convective
	667	lcc3dstra(i, k) = pclc(i, k)radocond(i, k)phase3d(i, k)
[1992]	668	lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
	669	lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
[3121]	670	!FC pour la glace (CAUSES)
[4704]	671	icc3dcon(i, k) = rnebcon(i, k)(1-phase3d(i, k))ccwcon(i, k) ! glace convective
	672	icc3dstra(i, k)= pclc(i, k)radocond(i, k)(1-phase3d(i, k))
[3121]	673	icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
	674	icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
	675	!FC (CAUSES)
	676
[1992]	677	! Compute cloud droplet radius as above in meter
[4704]	678	radius = 1.1((radocond(i,k)pplay(i,k)/(rdtemp(i,k)))/(4./3rpi1000. &
[1992]	679	cdnc(i,k)))**(1./3.)
	680	radius = max(radius, 5.E-6)
	681	! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
	682	reffclwc(i, k) = radius
	683	reffclwc(i, k) = reffclwc(i, k)*lcc3dcon(i, k)
	684	! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
	685	reffclws(i, k) = radius
	686	reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
[3274]	687	ENDDO !klev
	688	ENDDO !klon
[524]	689
[1992]	690	! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
[524]	691
[1992]	692	DO i = 1, klon
	693	cldnvi(i) = 0.
	694	lcc_integrat(i) = 0.
	695	height(i) = 0.
[1989]	696	DO k = 1, klev
[1992]	697	cldnvi(i) = cldnvi(i) + cdnc(i, k)lcc3d(i, k)dh(i, k)
	698	lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
	699	height(i) = height(i) + dh(i, k)
[3274]	700	ENDDO ! klev
[1992]	701	lcc_integrat(i) = lcc_integrat(i)/height(i)
	702	IF (lcc_integrat(i)<=1.0E-03) THEN
	703	cldnvi(i) = cldnvi(i)*lcc(i)/seuil_neb
	704	ELSE
	705	cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
[3274]	706	ENDIF
	707	ENDDO ! klon
[1337]	708
[1992]	709	DO i = 1, klon
	710	DO k = 1, klev
	711	IF (scdnc(i,k)<=0.0) scdnc(i, k) = 0.0
	712	IF (reffclws(i,k)<=0.0) reffclws(i, k) = 0.0
	713	IF (reffclwc(i,k)<=0.0) reffclwc(i, k) = 0.0
	714	IF (lcc3d(i,k)<=0.0) lcc3d(i, k) = 0.0
	715	IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
	716	IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
[3121]	717	!FC (CAUSES)
	718	IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
	719	IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
	720	!FC (CAUSES)
[3274]	721	ENDDO
[1992]	722	IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
	723	IF (cldncl(i)<=0.0) cldncl(i) = 0.0
	724	IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
	725	IF (lcc(i)<=0.0) lcc(i) = 0.0
[3274]	726	ENDDO
[1337]	727
[3274]	728	ENDIF !ok_cdnc
[1337]	729
[3245]	730	first=.false. !to be sure
	731
[1992]	732	RETURN
[1337]	733
[4683]	734	END SUBROUTINE cloud_optics_prop
	735
	736	END MODULE lmdz_cloud_optics_prop

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format