Context Navigation

MISR_simulator.F @ 4106

Last change on this file since 4106 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

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

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

Download in other formats:

Original Format