Context Navigation

icarus.F90 @ 5218

Last change on this file since 5218 was 5185, checked in by abarral, 3 months ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 30.1 KB

Line
1	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2	! Copyright (c) 2009, Lawrence Livemore National Security Limited Liability
3	! All rights reserved.
4
5	! Redistribution and use in source and binary forms, with or without modification, are
6	! permitted provided that the following conditions are met:
7
8	! 1. Redistributions of source code must retain the above copyright notice, this list of
9	! conditions and the following disclaimer.
10
11	! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12	! of conditions and the following disclaimer in the documentation and/or other
13	! materials provided with the distribution.
14
15	! 3. Neither the name of the copyright holder nor the names of its contributors may be
16	! used to endorse or promote products derived from this software without specific prior
17	! written permission.
18
19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20	! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21	! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22	! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24	! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26	! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28
29	! History
30	! May 2015 - D. Swales - Modified for COSPv2.0
31	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32	MODULE MOD_ICARUS
33	USE COSP_KINDS, ONLY: wp
34	USE COSP_PHYS_CONSTANTS, ONLY: amd,amw,avo,grav
35	use MOD_COSP_STATS, ONLY: hist2D
36	USE MOD_COSP_CONFIG, ONLY: R_UNDEF,numISCCPTauBins,numISCCPPresBins,isccp_histTau, &
37	isccp_histPres
38	implicit none
39
40	! Shared Parameters
41	integer,parameter :: &
42	ncolprint = 0 ! Flag for debug printing (set as parameter to increase performance)
43
44	! Cloud-top height determination
45	integer :: &
46	isccp_top_height, & ! Top height adjustment method
47	isccp_top_height_direction ! Direction for finding atmosphere pressure level
48
49	! Parameters used by icarus
50	real(wp),parameter :: &
51	tauchk = -1._wp*log(0.9999999_wp), & ! Lower limit on optical depth
52	isccp_taumin = 0.3_wp, & ! Minimum optical depth for joint-hostogram
53	pstd = 1013250._wp, & ! Mean sea-level pressure (Pa)
54	isccp_t0 = 296._wp, & ! Mean surface temperature (K)
55	output_missing_value = -1.E+30 ! Missing values
56
57	contains
58	! ##########################################################################
59	! ##########################################################################
60	SUBROUTINE ICARUS(debug,debugcol,npoints,sunlit,nlev,ncol,pfull, &
61	phalf,qv,cc,conv,dtau_s,dtau_c,th,thd,frac_out,skt,emsfc_lw,at,&
62	dem_s,dem_c,fq_isccp,totalcldarea, meanptop,meantaucld, &
63	meanalbedocld, meantb,meantbclr,boxtau,boxptop,levmatch)
64
65	! INPUTS
66	INTEGER,intent(in) :: & !
67	npoints, & ! Number of model points in the horizontal
68	nlev, & ! Number of model levels in column
69	ncol, & ! Number of subcolumns
70	debug, & ! Debug flag
71	debugcol ! Debug column flag
72	INTEGER,intent(in),dimension(npoints) :: & !
73	sunlit ! 1 for day points, 0 for night time
74	REAL(WP),intent(in) :: & !
75	emsfc_lw ! 10.5 micron emissivity of surface (fraction)
76	REAL(WP),intent(in),dimension(npoints) :: & !
77	skt ! Skin Temperature (K)
78	REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: & !
79	frac_out ! Boxes gridbox divided up into subcolumns
80	REAL(WP),intent(in),dimension(npoints,nlev) :: & !
81	pfull, & ! Pressure of full model levels (Pascals)
82	qv, & ! Water vapor specific humidity (kg vapor/ kg air)
83	cc, & ! Cloud cover in each model level (fraction)
84	conv, & ! Convective cloud cover in each model
85	at, & ! Temperature in each model level (K)
86	dem_c, & ! Emissivity for convective clouds
87	dem_s, & ! Emissivity for stratiform clouds
88	dtau_c, & ! Optical depth for convective clouds
89	dtau_s ! Optical depth for stratiform clouds
90	REAL(WP),intent(in),dimension(npoints,nlev+1) :: & !
91	phalf ! Pressure of half model levels (Pascals)!
92	integer,intent(in) :: th,thd
93
94	! OUTPUTS
95	REAL(WP),intent(out),dimension(npoints,7,7) :: &
96	fq_isccp ! The fraction of the model grid box covered by clouds
97	REAL(WP),intent(out),dimension(npoints) :: &
98	totalcldarea, & ! The fraction of model grid box columns with cloud present
99	meanptop, & ! Mean cloud top pressure (mb) - linear averaging
100	meantaucld, & ! Mean optical thickness
101	meanalbedocld, & ! Mean cloud albedo
102	meantb, & ! Mean all-sky 10.5 micron brightness temperature
103	meantbclr ! Mean clear-sky 10.5 micron brightness temperature
104	REAL(WP),intent(out),dimension(npoints,ncol) :: &
105	boxtau, & ! Optical thickness in each column
106	boxptop ! Cloud top pressure (mb) in each column
107	INTEGER,intent(out),dimension(npoints,ncol) :: &
108	levmatch ! Used for icarus unit testing only
109
110
111	! INTERNAL VARIABLES
112	CHARACTER(len=10) :: ftn09
113	REAL(WP),dimension(npoints,ncol) :: boxttop
114	REAL(WP),dimension(npoints,ncol,nlev) :: dtau,demIN
115	INTEGER :: j,ilev,ibox
116	INTEGER,dimension(nlev,ncol ) :: acc
117
118	! PARAMETERS
119	character ,parameter, dimension(6) :: cchar=(/' ','-','1','+','I','+'/)
120	character(len=1),parameter,dimension(6) :: cchar_realtops=(/ ' ',' ','1','1','I','I'/)
121	! ##########################################################################
122
123	call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dem_c,dem_s,demIN)
124	call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dtau_c,dtau_s,dtau)
125
126	call ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demIN,skt,emsfc_lw,qv,at, &
127	pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
128
129	call ICARUS_COLUMN(npoints,ncol,boxtau,boxptop/100._wp,sunlit,boxttop,&
130	fq_isccp,meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
131
132	! ##########################################################################
133	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
134	! ##########################################################################
135
136	if (debugcol.ne.0) then
137	DO j=1,npoints,debugcol
138
139	! Produce character output
140	DO ilev=1,nlev
141	acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2
142	where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
143	enddo
144
145	write(ftn09,11) j
146	11 format('ftn09.',i4.4)
147	open(9, FILE=ftn09, FORM='FORMATTED')
148
149	write(9,'(a1)') ' '
150	write(9,'(10i5)') (ilev,ilev=5,nlev,5)
151	write(9,'(a1)') ' '
152
153	DO ibox=1,ncol
154	write(9,'(40(a1),1x,40(a1))') &
155	(cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),&
156	(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
157	end do
158	close(9)
159
160	enddo
161	end if
162
163	return
164	END SUBROUTINE ICARUS
165
166	! ############################################################################
167	! ############################################################################
168	! ############################################################################
169	SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,at, &
170	pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
171	! Inputs
172	INTEGER, intent(in) :: &
173	ncol, & ! Number of subcolumns
174	npoints, & ! Number of horizontal gridpoints
175	nlev ! Number of vertical levels
176	INTEGER, intent(in), dimension(npoints) :: &
177	sunlit ! 1=day 0=night
178	REAL(WP),intent(in) :: &
179	emsfc_lw ! 10.5 micron emissivity of surface (fraction)
180	REAL(WP),intent(in), dimension(npoints) :: &
181	skt ! Skin temperature
182	REAL(WP),intent(in), dimension(npoints,nlev) :: &
183	at, & ! Temperature
184	pfull, & ! Presure
185	qv ! Specific humidity
186	REAL(WP),intent(in), dimension(npoints,ncol,nlev) :: &
187	frac_out, & ! Subcolumn cloud cover
188	dtau, & ! Subcolumn optical thickness
189	demIN ! Subcolumn emissivity
190	REAL(WP),intent(in), dimension(npoints,nlev+1) :: &
191	phalf ! Pressure at model half levels
192
193	! Outputs
194	REAL(WP),intent(inout),dimension(npoints) :: &
195	meantbclr ! Mean clear-sky 10.5 micron brightness temperature
196	REAL(WP),intent(inout),dimension(npoints,ncol) :: &
197	boxtau, & ! Optical thickness in each column
198	boxptop, & ! Cloud top pressure (mb) in each column
199	boxttop ! Cloud top temperature in each column
200	INTEGER, intent(inout),dimension(npoints,ncol) :: levmatch
201
202	! Local Variables
203	INTEGER :: &
204	j,ibox,ilev,k1,k2,icycle
205	INTEGER,dimension(npoints) :: &
206	nmatch,itrop
207	INTEGER,dimension(npoints,nlev-1) :: &
208	match
209	REAL(WP) :: &
210	logp,logp1,logp2,atd
211	REAL(WP),dimension(npoints) :: &
212	bb,attropmin,attrop,ptrop,atmax,btcmin,transmax,tauir,taumin,fluxtopinit,press, &
213	dpress,atmden,rvh20,rhoave,rh20s,rfrgn,tmpexp,tauwv,wk,trans_layers_above_clrsky, &
214	fluxtop_clrsky
215	REAL(WP),dimension(npoints,nlev) :: &
216	dem_wv
217	REAL(WP),dimension(npoints,ncol) :: &
218	trans_layers_above,dem,tb,emcld,fluxtop,tau,ptop
219
220	! ####################################################################################
221	! Compute cloud optical depth for each column by summing up subcolumns
222	tau(1:npoints,1:ncol) = 0._wp
223	tau(1:npoints,1:ncol) = sum(dtau,dim=3)
224
225	! Set tropopause values
226	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
227	ptrop(1:npoints) = 5000._wp
228	attropmin(1:npoints) = 400._wp
229	atmax(1:npoints) = 0._wp
230	attrop(1:npoints) = 120._wp
231	itrop(1:npoints) = 1
232
233	DO ilev=1,nlev
234	where(pfull(1:npoints,ilev) .lt. 40000. .AND. &
235	pfull(1:npoints,ilev) .gt. 5000. .AND. &
236	at(1:npoints,ilev) .lt. attropmin(1:npoints))
237	ptrop(1:npoints) = pfull(1:npoints,ilev)
238	attropmin(1:npoints) = at(1:npoints,ilev)
239	attrop(1:npoints) = attropmin(1:npoints)
240	itrop = ilev
241	endwhere
242	enddo
243
244	DO ilev=1,nlev
245	atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),&
246	at(1:npoints,ilev) .gt. atmax(1:npoints) .AND. ilev .ge. itrop(1:npoints))
247	enddo
248	end if
249
250	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
251	! ############################################################################
252	! Clear-sky radiance calculation
253
254	! Compute water vapor continuum emissivity this treatment follows Schwarkzopf
255	! and Ramasamy JGR 1999,vol 104, pages 9467-9499. The emissivity is calculated
256	! at a wavenumber of 955 cm-1, or 10.47 microns
257	! ############################################################################
258	DO ilev=1,nlev
259	press(1:npoints) = pfull(1:npoints,ilev)*10._wp
260	dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10
261	atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp)
262	rvh20(1:npoints) = qv(1:npoints,ilev)*amd/amw
263	wk(1:npoints) = rvh20(1:npoints)avoatmden/amd
264	rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev))
265	rh20s(1:npoints) = rvh20(1:npoints)*rhoave(1:npoints)
266	rfrgn(1:npoints) = rhoave(1:npoints)-rh20s(1:npoints)
267	tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0))
268	tauwv(1:npoints) = wk(1:npoints)1.e-20((0.0224697_wprh20s(1:npoints) &
269	tmpexp(1:npoints))+(3.41817e-7rfrgn(1:npoints)))0.98_wp
270	dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints))
271	enddo
272
273	fluxtop_clrsky(1:npoints) = 0._wp
274	trans_layers_above_clrsky(1:npoints) = 1._wp
275	DO ilev=1,nlev
276	! Black body emission at temperature of the layer
277	bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp )
278
279	! Increase TOA flux by flux emitted from layer times total transmittance in layers above
280	fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
281	dem_wv(1:npoints,ilev)bb(1:npoints)trans_layers_above_clrsky(1:npoints)
282
283	! Update trans_layers_above with transmissivity from this layer for next time around loop
284	trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*&
285	(1.-dem_wv(1:npoints,ilev))
286	enddo
287
288	! Add in surface emission
289	bb(1:npoints) = 1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
290	fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
291	emsfc_lw * bb(1:npoints)*trans_layers_above_clrsky(1:npoints)
292
293	! Clear Sky brightness temperature
294	meantbclr(1:npoints) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop_clrsky(1:npoints))))
295
296	! #################################################################################
297	! All-sky radiance calculation
298	! #################################################################################
299
300	fluxtop(1:npoints,1:ncol) = 0._wp
301	trans_layers_above(1:npoints,1:ncol) = 1._wp
302	DO ilev=1,nlev
303	! Black body emission at temperature of the layer
304	bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp)
305
306	DO ibox=1,ncol
307	! Emissivity
308	dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
309	1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
310	demIN(1:npoints,ibox,ilev) .eq. 0)
311
312	! Increase TOA flux emitted from layer
313	fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)bbtrans_layers_above(1:npoints,ibox)
314
315	! Update trans_layer by emitted layer from above
316	trans_layers_above(1:npoints,ibox) = trans_layers_above(1:npoints,ibox)*(1._wp-dem(1:npoints,ibox))
317	enddo
318	enddo
319
320	! Add in surface emission
321	bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
322	DO ibox=1,ncol
323	fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lwbb(1:npoints)trans_layers_above(1:npoints,ibox)
324	end do
325
326	! All Sky brightness temperature
327	boxttop(1:npoints,1:ncol) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop(1:npoints,1:ncol))))
328
329	! #################################################################################
330	! Cloud-Top Temperature
331
332	! Now that you have the top of atmosphere radiance, account for ISCCP
333	! procedures to determine cloud top temperature account for partially
334	! transmitting cloud recompute flux ISCCP would see assuming a single layer
335	! cloud. NOTE choice here of 2.13, as it is primarily ice clouds which have
336	! partial emissivity and need the adjustment performed in this section. If it
337	! turns out that the cloud brightness temperature is greater than 260K, then
338	! the liquid cloud conversion factor of 2.56 is used. NOTE that this is
339	! discussed on pages 85-87 of the ISCCP D level documentation
340	! (Rossow et al. 1996)
341	! #################################################################################
342
343	! Compute minimum brightness temperature and optical depth
344	btcmin(1:npoints) = 1._wp / ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp )
345
346	DO ibox=1,ncol
347	transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints))
348	tauir(1:npoints) = tau(1:npoints,ibox)/2.13_wp
349	taumin(1:npoints) = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
350	if (isccp_top_height .eq. 1) then
351	DO j=1,npoints
352	if (transmax(j) .gt. 0.001 .AND. transmax(j) .le. 0.9999999) then
353	fluxtopinit(j) = fluxtop(j,ibox)
354	tauir(j) = tau(j,ibox)/2.13_wp
355	endif
356	enddo
357	DO icycle=1,2
358	DO j=1,npoints
359	if (tau(j,ibox) .gt. (tauchk)) then
360	if (transmax(j) .gt. 0.001 .AND. transmax(j) .le. 0.9999999) then
361	emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j) )
362	fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
363	fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
364	tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
365	if (tb(j,ibox) .gt. 260.) then
366	tauir(j) = tau(j,ibox) / 2.56_wp
367	end if
368	end if
369	end if
370	enddo
371	enddo
372	endif
373
374	! Cloud-top temperature
375	where(tau(1:npoints,ibox) .gt. tauchk)
376	tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox))))
377	where (isccp_top_height .eq. 1 .AND. tauir(1:npoints) .lt. taumin(1:npoints))
378	tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp
379	tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
380	endwhere
381	endwhere
382
383	! Clear-sky brightness temperature
384	where(tau(1:npoints,ibox) .le. tauchk)
385	tb(1:npoints,ibox) = meantbclr(1:npoints)
386	endwhere
387	enddo
388	else
389	meantbclr(1:npoints) = output_missing_value
390	end if
391
392	! ####################################################################################
393	! Cloud-Top Pressure
394
395	! The 2 methods differ according to whether or not you use the physical cloud
396	! top pressure (isccp_top_height = 2) or the radiatively determined cloud top
397	! pressure (isccp_top_height = 1 or 3)
398	! ####################################################################################
399	DO ibox=1,ncol
400	!segregate according to optical thickness
401	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
402
403	! Find level whose temperature most closely matches brightness temperature
404	nmatch(1:npoints)=0
405	DO k1=1,nlev-1
406	ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)
407	DO j=1,npoints
408	if (ilev .ge. itrop(j) .AND. &
409	((at(j,ilev) .ge. tb(j,ibox) .AND. &
410	at(j,ilev+1) .le. tb(j,ibox)) .or. &
411	(at(j,ilev) .le. tb(j,ibox) .AND. &
412	at(j,ilev+1) .ge. tb(j,ibox)))) then
413	nmatch(j)=nmatch(j)+1
414	match(j,nmatch(j))=ilev
415	endif
416	enddo
417	enddo
418
419	DO j=1,npoints
420	if (nmatch(j) .ge. 1) then
421	k1 = match(j,nmatch(j))
422	k2 = k1 + 1
423	logp1 = log(pfull(j,k1))
424	logp2 = log(pfull(j,k2))
425	atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
426	logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
427	ptop(j,ibox) = exp(logp)
428	levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox)))
429	else
430	if (tb(j,ibox) .le. attrop(j)) then
431	ptop(j,ibox)=ptrop(j)
432	levmatch(j,ibox)=itrop(j)
433	end if
434	if (tb(j,ibox) .ge. atmax(j)) then
435	ptop(j,ibox)=pfull(j,nlev)
436	levmatch(j,ibox)=nlev
437	end if
438	end if
439	enddo
440	else
441	ptop(1:npoints,ibox)=0.
442	DO ilev=1,nlev
443	where((ptop(1:npoints,ibox) .eq. 0. ) .AND.(frac_out(1:npoints,ibox,ilev) .ne. 0))
444	ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
445	levmatch(1:npoints,ibox)=ilev
446	endwhere
447	end do
448	end if
449	where(tau(1:npoints,ibox) .le. tauchk)
450	ptop(1:npoints,ibox)=0._wp
451	levmatch(1:npoints,ibox)=0._wp
452	endwhere
453	enddo
454
455	! ####################################################################################
456	! Compute subcolumn pressure and optical depth
457	! ####################################################################################
458	boxtau(1:npoints,1:ncol) = output_missing_value
459	boxptop(1:npoints,1:ncol) = output_missing_value
460	DO ibox=1,ncol
461	DO j=1,npoints
462	if (tau(j,ibox) .gt. (tauchk) .AND. ptop(j,ibox) .gt. 0.) then
463	if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
464	boxtau(j,ibox) = tau(j,ibox)
465	boxptop(j,ibox) = ptop(j,ibox)!/100._wp
466	endif
467	endif
468	enddo
469	enddo
470
471	END SUBROUTINE ICARUS_SUBCOLUMN
472
473	! ######################################################################################
474	! SUBROUTINE icarus_column
475	! ######################################################################################
476	SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp, &
477	meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
478	! Inputs
479	INTEGER, intent(in) :: &
480	ncol, & ! Number of subcolumns
481	npoints ! Number of horizontal gridpoints
482	INTEGER, intent(in),dimension(npoints) :: &
483	sunlit ! day=1 night=0
484	REAL(WP),intent(in),dimension(npoints,ncol) :: &
485	boxttop, & ! Subcolumn top temperature
486	boxptop, & ! Subcolumn cloud top pressure
487	boxtau ! Subcolumn optical depth
488
489	! Outputs
490	REAL(WP),intent(inout),dimension(npoints) :: &
491	meanalbedocld, & ! Gridmean cloud albedo
492	meanptop, & ! Gridmean cloud top pressure (mb) - linear averaging
493	meantaucld, & ! Gridmean optical thickness
494	totalcldarea, & ! The fraction of model grid box columns with cloud present
495	meantb ! Gridmean all-sky 10.5 micron brightness temperature
496	REAL(WP),intent(inout),dimension(npoints,7,7) :: &
497	fq_isccp ! The fraction of the model grid box covered by clouds
498
499	! Local Variables
500	INTEGER :: j,ilev,ilev2
501	REAL(WP),dimension(npoints,ncol) :: albedocld
502	LOGICAL, dimension(npoints,ncol) :: box_cloudy
503
504	! Variables for new joint-histogram implementation
505	logical,dimension(ncol) :: box_cloudy2
506
507	! ####################################################################################
508	! Brightness Temperature
509	! ####################################################################################
510	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
511	meantb(1:npoints)=sum(boxttop,2)/ncol
512	else
513	meantb(1:npoints) = output_missing_value
514	endif
515
516	! ####################################################################################
517	! Determines ISCCP cloud type frequencies
518
519	! Now that boxptop and boxtau have been determined, determine amount of each of the
520	! 49 ISCCP cloud types. Also compute grid box mean cloud top pressure and
521	! optical thickness. The mean cloud top pressure and optical thickness are
522	! averages over the cloudy area only. The mean cloud top pressure is a linear
523	! average of the cloud top pressures. The mean cloud optical thickness is
524	! computed by converting optical thickness to an albedo, averaging in albedo
525	! units, then converting the average albedo back to a mean optical thickness.
526	! ####################################################################################
527
528	! Initialize
529	albedocld(1:npoints,1:ncol) = 0._wp
530	box_cloudy(1:npoints,1:ncol) = .false.
531
532	! Reset frequencies
533	!fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7)
534	DO ilev=1,7
535	DO ilev2=1,7
536	DO j=1,npoints !
537	if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
538	fq_isccp(j,ilev,ilev2)= 0.
539	else
540	fq_isccp(j,ilev,ilev2)= output_missing_value
541	end if
542	enddo
543	enddo
544	enddo
545
546
547	! Reset variables need for averaging cloud properties
548	where(sunlit .eq. 1 .or. isccp_top_height .eq. 3)
549	totalcldarea(1:npoints) = 0._wp
550	meanalbedocld(1:npoints) = 0._wp
551	meanptop(1:npoints) = 0._wp
552	meantaucld(1:npoints) = 0._wp
553	elsewhere
554	totalcldarea(1:npoints) = output_missing_value
555	meanalbedocld(1:npoints) = output_missing_value
556	meanptop(1:npoints) = output_missing_value
557	meantaucld(1:npoints) = output_missing_value
558	endwhere
559
560	! Compute column quantities and joint-histogram
561	DO j=1,npoints
562	! Subcolumns that are cloudy(true) and not(false)
563	box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .AND. boxptop(j,1:ncol) .gt. 0.)
564
565	! Compute joint histogram and column quantities for points that are sunlit and cloudy
566	if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then
567	! Joint-histogram
568	call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, &
569	isccp_histPres,numISCCPPresBins,fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins))
570	fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins) = &
571	fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol
572
573	! Column cloud area
574	totalcldarea(j) = real(count(box_cloudy2(1:ncol) .AND. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol
575
576	! Subcolumn cloud albedo
577	!albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)0.895_wp)/((boxtau(j,1:ncol)0.895_wp)+6.82_wp),&
578	! 0._wp,box_cloudy2(1:ncol) .AND. boxtau(j,1:ncol) .gt. isccp_taumin)
579	where(box_cloudy2(1:ncol) .AND. boxtau(j,1:ncol) .gt. isccp_taumin)
580	albedocld(j,1:ncol) = (boxtau(j,1:ncol)0.895_wp)/((boxtau(j,1:ncol)0.895_wp)+6.82_wp)
581	elsewhere
582	albedocld(j,1:ncol) = 0._wp
583	endwhere
584
585	! Column cloud albedo
586	meanalbedocld(j) = sum(albedocld(j,1:ncol))/ncol
587
588	! Column cloud top pressure
589	meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .AND. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol
590	endif
591	enddo
592
593	! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
594	where(totalcldarea(1:npoints) .gt. 0)
595	meanptop(1:npoints) = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
596	meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
597	meantaucld(1:npoints) = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp)
598	elsewhere
599	meanptop(1:nPoints) = output_missing_value
600	meanalbedocld(1:nPoints) = output_missing_value
601	meantaucld(1:nPoints) = output_missing_value
602	endwhere
603	!meanptop(1:npoints) = merge(100._wp*meanptop(1:npoints)/totalcldarea(1:npoints),&
604	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
605	!meanalbedocld(1:npoints) = merge(meanalbedocld(1:npoints)/totalcldarea(1:npoints), &
606	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
607	!meantaucld(1:npoints) = merge((6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp), &
608	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
609
610	! Represent in percent
611	where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp
612	where(fq_isccp .ne. output_missing_value) fq_isccp = fq_isccp*100._wp
613
614
615	END SUBROUTINE ICARUS_column
616
617	subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
618	! INPUTS
619	integer,intent(in) :: &
620	dim1, & ! Dimension 1 extent (Horizontal)
621	dim2, & ! Dimension 2 extent (Subcolumn)
622	dim3 ! Dimension 3 extent (Vertical)
623	real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
624	flag ! Logical to determine the of merge var1IN and var2IN
625	real(wp),intent(in),dimension(dim1, dim3) :: &
626	varIN1, & ! Input field 1
627	varIN2 ! Input field 2
628	! OUTPUTS
629	real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
630	varOUT ! Merged output field
631	! LOCAL VARIABLES
632	integer :: j
633
634	varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
635	DO j=1,dim2
636	where(flag(:,j,:) .eq. 1)
637	varOUT(:,j,:) = varIN2
638	endwhere
639	where(flag(:,j,:) .eq. 2)
640	varOUT(:,j,:) = varIN1
641	endwhere
642	enddo
643	end subroutine cosp_simulator_optics
644	end module MOD_ICARUS
645

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

Download in other formats:

Original Format