Context Navigation

MISR_simulator.F @ 5136

Last change on this file since 5136 was 5099, checked in by abarral, 4 months ago
Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90 Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad Create lmdz_yoerad.f90 (lint) Remove "!" on otherwise empty line
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
File size: 16.0 KB

Line
1
2	! Copyright (c) 2009, Roger Marchand, version 1.2
3	! All rights reserved.
4	! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
5	! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MISR_simulator/MISR_simulator.f $
6
7	! Redistribution and use in source and binary forms, with or without modification, are permitted
8	! provided that the following conditions are met:
9
10	! * Redistributions of source code must retain the above copyright notice, this list of
11	! conditions and the following disclaimer.
12	! * Redistributions in binary form must reproduce the above copyright notice, this list
13	! of conditions and the following disclaimer in the documentation and/or other materials
14	! provided with the distribution.
15	! * Neither the name of the University of Washington nor the names of its contributors may be used
16	! to endorse or promote products derived from this software without specific prior written permission.
17
18	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
19	! BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20	! SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23	! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25	SUBROUTINE MISR_simulator(
26	& npoints,
27	& nlev,
28	& ncol,
29	& sunlit,
30	& zfull,
31	& at,
32	& dtau_s,
33	& dtau_c,
34	& frac_out,
35	& missing_value,
36	& fq_MISR_TAU_v_CTH,
37	& dist_model_layertops,
38	& MISR_mean_ztop,
39	& MISR_cldarea
40	& )
41
42
43	implicit none
44	integer n_MISR_CTH
45	parameter(n_MISR_CTH=16)
46
47	! -----
48	! Input
49	! -----
50
51	INTEGER npoints ! if ncol ==1, the number of model points in the horizontal grid
52	! else the number of GCM grid points
53
54	INTEGER nlev ! number of model vertical levels
55
56	INTEGER ncol ! number of model sub columns
57	! (must already be generated in via scops and passed to this
58	! routine via the variable frac_out )
59
60	INTEGER sunlit(npoints) ! 1 for day points, 0 for night time
61
62	REAL zfull(npoints,nlev) ! height (in meters) of full model levels (i.e. midpoints)
63	! zfull(npoints,1) is top level of model
64	! zfull(npoints,nlev) is bottom level of model (closest point to surface)
65
66	REAL at(npoints,nlev) ! temperature in each model level (K)
67
68	REAL dtau_s(npoints,nlev) ! visible wavelength cloud optical depth ... for "stratiform" condensate
69	! NOTE: this the cloud optical depth of only the
70	! the model cell (i,j)
71
72	REAL dtau_c(npoints,nlev) ! visible wavelength cloud optical depth ... for "convective" condensate
73	! NOTE: this the cloud optical depth of only the
74	! the model cell (i,j)
75
76	REAL frac_out(npoints,ncol,nlev) ! NOTE: only need if columns>1 ... subgrid scheme in use.
77
78	REAL missing_value
79
80	! ------
81	! Outputs
82	! ------
83
84	REAL fq_MISR_TAU_v_CTH(npoints,7,n_MISR_CTH)
85	REAL dist_model_layertops(npoints,n_MISR_CTH)
86	REAL MISR_cldarea(npoints) ! fractional area coverged by clouds
87	REAL MISR_mean_ztop(npoints) ! mean cloud top hieght(m) MISR would observe
88	! NOTE: == 0 if area ==0
89
90
91	! ------
92	! Working variables
93	! ------
94
95	REAL tau(npoints,ncol) ! total column optical depth ...
96
97	INTEGER j,ilev,ilev2,ibox,k
98	INTEGER itau
99
100	LOGICAL box_cloudy(npoints,ncol)
101
102	real isccp_taumin
103	real boxarea
104	real tauchk
105	REAL box_MISR_ztop(npoints,ncol) ! cloud top hieght(m) MISR would observe
106
107	integer thres_crossed_MISR
108	integer loop,iMISR_ztop
109
110	real dtau, cloud_dtau, MISR_penetration_height,ztest
111
112	real MISR_CTH_boundaries(n_MISR_CTH+1)
113
114	DATA MISR_CTH_boundaries / -99, 0, 0.5, 1, 1.5, 2, 2.5, 3,
115	c 4, 5, 7, 9, 11, 13, 15, 17, 99 /
116
117	DATA isccp_taumin / 0.3 /
118
119	tauchk = -1.*log(0.9999999)
120
121	! For each GCM cell or horizontal model grid point ...
122
123	do j=1,npoints
124
125	! estimate distribution of Model layer tops
126
127	dist_model_layertops(j,:)=0
128
129	do ilev=1,nlev
130
131	! define location of "layer top"
132	if(ilev.eq.1 .or. ilev.eq.nlev) then
133	ztest=zfull(j,ilev)
134	else
135	ztest=0.5*(zfull(j,ilev)+zfull(j,ilev-1))
136	endif
137
138	! find MISR layer that contains this level
139	! note, the first MISR level is "no height" level
140	iMISR_ztop=2
141	do loop=2,n_MISR_CTH
142
143	if ( ztest .gt.
144	& 1000*MISR_CTH_boundaries(loop+1) ) then
145
146	iMISR_ztop=loop+1
147	endif
148	enddo
149
150	dist_model_layertops(j,iMISR_ztop)=
151	& dist_model_layertops(j,iMISR_ztop)+1
152	enddo
153
154	! compute total cloud optical depth for each column
155
156	do ibox=1,ncol
157
158	! Initialize tau to zero in each subcolum
159	tau(j,ibox)=0.
160	box_cloudy(j,ibox)=.false.
161	box_MISR_ztop(j,ibox)=0
162
163	! initialize threshold detection for each sub column
164	thres_crossed_MISR=0;
165
166	do ilev=1,nlev
167
168	dtau=0
169
170	if (frac_out(j,ibox,ilev).eq.1) then
171	dtau = dtau_s(j,ilev)
172	endif
173
174	if (frac_out(j,ibox,ilev).eq.2) then
175	dtau = dtau_c(j,ilev)
176	end if
177
178	tau(j,ibox)=tau(j,ibox)+ dtau
179
180
181	! NOW for MISR ..
182	! if there a cloud ... start the counter ... store this height
183	if(thres_crossed_MISR .eq. 0 .and. dtau .gt. 0.) then
184
185	! first encountered a "cloud"
186	thres_crossed_MISR=1
187	cloud_dtau=0
188	endif
189
190	if( thres_crossed_MISR .lt. 99 .and.
191	& thres_crossed_MISR .gt. 0 ) then
192
193	if( dtau .eq. 0.) then
194
195	! we have come to the end of the current cloud
196	! layer without yet selecting a CTH boundary.
197	! ... restart cloud tau counter
198	cloud_dtau=0
199	else
200	! add current optical depth to count for
201	! the current cloud layer
202	cloud_dtau=cloud_dtau+dtau
203	endif
204
205	! if the cloud is continuous but optically thin (< 1)
206	! from above the current layer cloud top to the current level
207	! then MISR will like see a top below the top of the current
208	! layer
209	if( dtau.gt.0 .and. (cloud_dtau-dtau) .lt. 1) then
210
211	if(dtau .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then
212
213	! MISR will likely penetrate to some point
214	! within this layer ... the middle
215	MISR_penetration_height=zfull(j,ilev)
216
217	else
218	! take the OD = 1.0 level into this layer
219	MISR_penetration_height=
220	& 0.5*(zfull(j,ilev)+zfull(j,ilev-1)) -
221	& 0.5*(zfull(j,ilev-1)-zfull(j,ilev+1))
222	& /dtau
223	endif
224
225	box_MISR_ztop(j,ibox)=MISR_penetration_height
226
227	endif
228
229	! check for a distinctive water layer
230	if(dtau .gt. 1 .and. at(j,ilev).gt.273 ) then
231
232	! must be a water cloud ...
233	! take this as CTH level
234	thres_crossed_MISR=99
235	endif
236
237	! if the total column optical depth is "large" than
238	! MISR can't seen anything else ... set current point as CTH level
239	if(tau(j,ibox) .gt. 5) then
240
241	thres_crossed_MISR=99
242	endif
243
244	endif ! MISR CTH booundary not set
245
246	enddo !ilev - loop over vertical levesl
247
248	! written by roj 5/2006
249	! check to see if there was a cloud for which we didn't
250	! set a MISR cloud top boundary
251	if( thres_crossed_MISR .eq. 1) then
252
253	! if the cloud has a total optical depth of greater
254	! than ~ 0.5 MISR will still likely pick up this cloud
255	! with a height near the true cloud top
256	! otherwise there should be no CTH
257	if( tau(j,ibox) .gt. 0.5) then
258
259	! keep MISR detected CTH
260
261	elseif(tau(j,ibox) .gt. 0.2) then
262
263	! MISR may detect but wont likley have a good height
264	box_MISR_ztop(j,ibox)=-1
265
266	else
267	! MISR not likely to even detect.
268	! so set as not cloudy
269	box_MISR_ztop(j,ibox)=0
270
271	endif
272
273	endif
274
275	enddo ! loop of subcolumns
276	enddo ! loop of gridpoints
277
278	! Modify MISR CTH for satellite spatial / pattern matcher effects
279
280	! Code in this region added by roj 5/2006 to account
281	! for spatial effect of the MISR pattern matcher.
282	! Basically, if a column is found between two neighbors
283	! at the same CTH, and that column has no hieght or
284	! a lower CTH, THEN misr will tend to but place the
285	! odd column at the same height as it neighbors.
286
287	! This setup assumes the columns represent a about a 1 to 4 km scale
288	! it will need to be modified significantly, otherwise
289	if(ncol.eq.1) then
290
291	! adjust based on neightboring points ... i.e. only 2D grid was input
292	do j=2,npoints-1
293
294	if(box_MISR_ztop(j-1,1).gt.0 .and.
295	& box_MISR_ztop(j+1,1).gt.0 ) then
296
297	if( abs( box_MISR_ztop(j-1,1) -
298	& box_MISR_ztop(j+1,1) ) .lt. 500
299	& .and.
300	& box_MISR_ztop(j,1) .lt.
301	& box_MISR_ztop(j+1,1) ) then
302
303	box_MISR_ztop(j,1) =
304	& box_MISR_ztop(j+1,1)
305	endif
306
307	endif
308	enddo
309	else
310
311	! adjust based on neighboring subcolumns ....
312	do ibox=2,ncol-1
313
314	if(box_MISR_ztop(1,ibox-1).gt.0 .and.
315	& box_MISR_ztop(1,ibox+1).gt.0 ) then
316
317	if( abs( box_MISR_ztop(1,ibox-1) -
318	& box_MISR_ztop(1,ibox+1) ) .lt. 500
319	& .and.
320	& box_MISR_ztop(1,ibox) .lt.
321	& box_MISR_ztop(1,ibox+1) ) then
322
323	box_MISR_ztop(1,ibox) =
324	& box_MISR_ztop(1,ibox+1)
325	endif
326
327	endif
328	enddo
329
330	endif
331
332	! DETERMINE CLOUD TYPE FREQUENCIES
333
334	! Now that ztop and tau have been determined,
335	! determine amount of each cloud type
336	boxarea=1./real(ncol)
337	do j=1,npoints
338
339	! reset frequencies -- modified loop structure, roj 5/2006
340	do ilev=1,7 ! "tau loop"
341	do ilev2=1,n_MISR_CTH
342	fq_MISR_TAU_v_CTH(j,ilev,ilev2)=0.
343	enddo
344	enddo
345
346	MISR_cldarea(j)=0.
347	MISR_mean_ztop(j)=0.
348
349	do ibox=1,ncol
350
351	if (tau(j,ibox) .gt. (tauchk)) then
352	box_cloudy(j,ibox)=.true.
353	endif
354
355	itau = 0
356
357	if (box_cloudy(j,ibox)) then
358
359	!determine optical depth category
360	if (tau(j,ibox) .lt. isccp_taumin) then
361	itau=1
362	else if (tau(j,ibox) .ge. isccp_taumin
363	& .and. tau(j,ibox) .lt. 1.3) then
364	itau=2
365	else if (tau(j,ibox) .ge. 1.3
366	& .and. tau(j,ibox) .lt. 3.6) then
367	itau=3
368	else if (tau(j,ibox) .ge. 3.6
369	& .and. tau(j,ibox) .lt. 9.4) then
370	itau=4
371	else if (tau(j,ibox) .ge. 9.4
372	& .and. tau(j,ibox) .lt. 23.) then
373	itau=5
374	else if (tau(j,ibox) .ge. 23.
375	& .and. tau(j,ibox) .lt. 60.) then
376	itau=6
377	else if (tau(j,ibox) .ge. 60.) then
378	itau=7
379	endif
380
381	endif
382
383	! update MISR histograms and summary metrics - roj 5/2005
384	if (sunlit(j).eq.1) then
385
386	!if cloudy added by roj 5/2005
387	if( box_MISR_ztop(j,ibox).eq.0) then
388
389	! no cloud detected
390	iMISR_ztop=0
391
392	elseif( box_MISR_ztop(j,ibox).eq.-1) then
393
394	! cloud can be detected but too thin to get CTH
395	iMISR_ztop=1
396
397	fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)=
398	& fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea
399
400	else
401
402	! determine index for MISR bin set
403
404	iMISR_ztop=2
405
406	do loop=2,n_MISR_CTH
407
408	if ( box_MISR_ztop(j,ibox) .gt.
409	& 1000*MISR_CTH_boundaries(loop+1) ) then
410
411	iMISR_ztop=loop+1
412
413	endif
414	enddo
415
416	if(box_cloudy(j,ibox)) then
417
418	! there is an isccp clouds so itau(j) is defined
419	fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)=
420	& fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea
421
422	else
423	! MISR CTH resolution is trying to fill in a
424	! broken cloud scene where there is no condensate.
425	! The MISR CTH-1D-OD product will only put in a cloud
426	! if the MISR cloud mask indicates cloud.
427	! therefore we will not include this column in the histogram
428	! in reality aerosoal and 3D effects or bright surfaces
429	! could fool the MISR cloud mask
430
431	! the alternative is to count as very thin cloud ??
432	! fq_MISR_TAU_v_CTH(1,iMISR_ztop)=
433	! & fq_MISR_TAU_v_CTH(1,iMISR_ztop) + boxarea
434	endif
435
436
437	MISR_mean_ztop(j)=MISR_mean_ztop(j)+
438	& box_MISR_ztop(j,ibox)*boxarea
439
440	MISR_cldarea(j)=MISR_cldarea(j) + boxarea
441
442	endif
443	else
444	! Set to issing data. A. Bodas - 14/05/2010
445	do loop=1,n_MISR_CTH
446	do k=1,7
447	fq_MISR_TAU_v_CTH(j,k,loop) = missing_value
448	enddo
449	dist_model_layertops(j,loop) = missing_value
450	enddo
451	MISR_cldarea(j) = missing_value
452	MISR_mean_ztop(npoints) = missing_value
453
454	endif ! is sunlight ?
455
456	enddo ! ibox - loop over subcolumns
457
458	if( MISR_cldarea(j) .gt. 0.) then
459	MISR_mean_ztop(j)= MISR_mean_ztop(j) / MISR_cldarea(j) ! roj 5/2006
460	endif
461
462	enddo ! loop over grid points
463
464	return
465	end

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

Download in other formats:

Original Format