My Project
 All Classes Files Functions Variables Macros
MISR_simulator.F
Go to the documentation of this file.
1 !
2 ! Copyright (c) 2009, Roger Marchand, version 1.2
3 ! All rights reserved.
4 !
5 ! Redistribution and use in source and binary forms, with or without modification, are permitted
6 ! provided that the following conditions are met:
7 !
8 ! * Redistributions of source code must retain the above copyright notice, this list of
9 ! conditions and the following disclaimer.
10 ! * Redistributions in binary form must reproduce the above copyright notice, this list
11 ! of conditions and the following disclaimer in the documentation and/or other materials
12 ! provided with the distribution.
13 ! * Neither the name of the University of Washington nor the names of its contributors may be used
14 ! to endorse or promote products derived from this software without specific prior written permission.
15 !
16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
17 ! BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
18 ! SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
20 ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
21 ! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
22 !
23 
24  SUBROUTINE misr_simulator(
25  & npoints,
26  & nlev,
27  & ncol,
28  & sunlit,
29  & zfull,
30  & at,
31  & dtau_s,
32  & dtau_c,
33  & frac_out,
34  & fq_misr_tau_v_cth,
35  & dist_model_layertops,
36  & misr_mean_ztop,
37  & misr_cldarea
38  & )
39 
40 
41  implicit none
42  integer n_misr_cth
43  parameter(n_misr_cth=16)
44 
45 ! -----
46 ! Input
47 ! -----
48 
49  INTEGER npoints ! if ncol ==1, the number of model points in the horizontal grid
50  ! else the number of GCM grid points
51 
52  INTEGER nlev ! number of model vertical levels
53 
54  INTEGER ncol ! number of model sub columns
55  ! (must already be generated in via scops and passed to this
56  ! routine via the variable frac_out )
57 
58  INTEGER sunlit(npoints) ! 1 for day points, 0 for night time
59 
60  REAL zfull(npoints,nlev) ! height (in meters) of full model levels (i.e. midpoints)
61  ! zfull(npoints,1) is top level of model
62  ! zfull(npoints,nlev) is bottom level of model (closest point to surface)
63 
64  REAL at(npoints,nlev) ! temperature in each model level (K)
65 
66  REAL dtau_s(npoints,nlev) ! visible wavelength cloud optical depth ... for "stratiform" condensate
67  ! NOTE: this the cloud optical depth of only the
68  ! the model cell (i,j)
69 
70  REAL dtau_c(npoints,nlev) ! visible wavelength cloud optical depth ... for "convective" condensate
71  ! NOTE: this the cloud optical depth of only the
72  ! the model cell (i,j)
73 
74  REAL frac_out(npoints,ncol,nlev) ! NOTE: only need if columns>1 ... subgrid scheme in use.
75 
76 ! ------
77 ! Outputs
78 ! ------
79 
80  REAL fq_misr_tau_v_cth(npoints,7,n_misr_cth)
81  REAL dist_model_layertops(npoints,n_misr_cth)
82  REAL misr_cldarea(npoints) ! fractional area coverged by clouds
83  REAL misr_mean_ztop(npoints) ! mean cloud top hieght(m) MISR would observe
84  ! NOTE: == 0 if area ==0
85 
86 
87 ! ------
88 ! Working variables
89 ! ------
90 
91  REAL tau(npoints,ncol) ! total column optical depth ...
92 
93  INTEGER j,ilev,ilev2,ibox
94  INTEGER itau
95 
96  LOGICAL box_cloudy(npoints,ncol)
97 
98  real isccp_taumin
99  real boxarea
100  real tauchk
101  REAL box_misr_ztop(npoints,ncol) ! cloud top hieght(m) MISR would observe
102 
103  integer thres_crossed_misr
104  integer loop,imisr_ztop
105 
106  real dtau, cloud_dtau, misr_penetration_height,ztest
107 
108  real misr_cth_boundaries(n_misr_cth+1)
109 
110  DATA misr_cth_boundaries / -99, 0, 0.5, 1, 1.5, 2, 2.5, 3,
111  c 4, 5, 7, 9, 11, 13, 15, 17, 99 /
112 
113  DATA isccp_taumin / 0.3 /
114 
115  tauchk = -1.*log(0.9999999)
116 
117  !
118  ! For each GCM cell or horizontal model grid point ...
119  !
120  do j=1,npoints
121 
122  !
123  ! estimate distribution of Model layer tops
124  !
125  dist_model_layertops(j,:)=0
126 
127  do ilev=1,nlev
128 
129  ! define location of "layer top"
130  if(ilev.eq.1 .or. ilev.eq.nlev) then
131  ztest=zfull(j,ilev)
132  else
133  ztest=0.5*(zfull(j,ilev)+zfull(j,ilev-1))
134  endif
135 
136  ! find MISR layer that contains this level
137  ! note, the first MISR level is "no height" level
138  imisr_ztop=2
139  do loop=2,n_misr_cth
140 
141  if ( ztest .gt.
142  & 1000*misr_cth_boundaries(loop+1) ) then
143 
144  imisr_ztop=loop+1
145  endif
146  enddo
147 
148  dist_model_layertops(j,imisr_ztop)=
149  & dist_model_layertops(j,imisr_ztop)+1
150  enddo
151 
152 
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 
279  !
280  ! Modify MISR CTH for satellite spatial / pattern matcher effects
281  !
282  ! Code in this region added by roj 5/2006 to account
283  ! for spatial effect of the MISR pattern matcher.
284  ! Basically, if a column is found between two neighbors
285  ! at the same CTH, and that column has no hieght or
286  ! a lower CTH, THEN misr will tend to but place the
287  ! odd column at the same height as it neighbors.
288  !
289  ! This setup assumes the columns represent a about a 1 to 4 km scale
290  ! it will need to be modified significantly, otherwise
291  if(ncol.eq.1) then
292 
293  ! adjust based on neightboring points ... i.e. only 2D grid was input
294  do j=2,npoints-1
295 
296  if(box_misr_ztop(j-1,1).gt.0 .and.
297  & box_misr_ztop(j+1,1).gt.0 ) then
298 
299  if( abs( box_misr_ztop(j-1,1) -
300  & box_misr_ztop(j+1,1) ) .lt. 500
301  & .and.
302  & box_misr_ztop(j,1) .lt.
303  & box_misr_ztop(j+1,1) ) then
304 
305  box_misr_ztop(j,1) =
306  & box_misr_ztop(j+1,1)
307  endif
308 
309  endif
310  enddo
311  else
312 
313  ! adjust based on neighboring subcolumns ....
314  do ibox=2,ncol-1
315 
316  if(box_misr_ztop(1,ibox-1).gt.0 .and.
317  & box_misr_ztop(1,ibox+1).gt.0 ) then
318 
319  if( abs( box_misr_ztop(1,ibox-1) -
320  & box_misr_ztop(1,ibox+1) ) .lt. 500
321  & .and.
322  & box_misr_ztop(1,ibox) .lt.
323  & box_misr_ztop(1,ibox+1) ) then
324 
325  box_misr_ztop(1,ibox) =
326  & box_misr_ztop(1,ibox+1)
327  endif
328 
329  endif
330  enddo
331 
332  endif
333 
334  !
335  ! DETERMINE CLOUD TYPE FREQUENCIES
336  !
337  ! Now that ztop and tau have been determined,
338  ! determine amount of each cloud type
339  boxarea=1./real(ncol)
340  do j=1,npoints
341 
342  ! reset frequencies -- modified loop structure, roj 5/2006
343  do ilev=1,7 ! "tau loop"
344  do ilev2=1,n_misr_cth
345  fq_misr_tau_v_cth(j,ilev,ilev2)=0.
346  enddo
347  enddo
348 
349  misr_cldarea(j)=0.
350  misr_mean_ztop(j)=0.
351 
352  do ibox=1,ncol
353 
354  if (tau(j,ibox) .gt. (tauchk)) then
355  box_cloudy(j,ibox)=.true.
356  endif
357 
358  itau = 0
359 
360  if (box_cloudy(j,ibox)) then
361 
362  !determine optical depth category
363  if (tau(j,ibox) .lt. isccp_taumin) then
364  itau=1
365  else if (tau(j,ibox) .ge. isccp_taumin
366  & .and. tau(j,ibox) .lt. 1.3) then
367  itau=2
368  else if (tau(j,ibox) .ge. 1.3
369  & .and. tau(j,ibox) .lt. 3.6) then
370  itau=3
371  else if (tau(j,ibox) .ge. 3.6
372  & .and. tau(j,ibox) .lt. 9.4) then
373  itau=4
374  else if (tau(j,ibox) .ge. 9.4
375  & .and. tau(j,ibox) .lt. 23.) then
376  itau=5
377  else if (tau(j,ibox) .ge. 23.
378  & .and. tau(j,ibox) .lt. 60.) then
379  itau=6
380  else if (tau(j,ibox) .ge. 60.) then
381  itau=7
382  endif
383 
384  endif
385 
386  ! update MISR histograms and summary metrics - roj 5/2005
387  if (sunlit(j).eq.1) then
388 
389  !if cloudy added by roj 5/2005
390  if( box_misr_ztop(j,ibox).eq.0) then
391 
392  ! no cloud detected
393  imisr_ztop=0
394 
395  elseif( box_misr_ztop(j,ibox).eq.-1) then
396 
397  ! cloud can be detected but too thin to get CTH
398  imisr_ztop=1
399 
400  fq_misr_tau_v_cth(j,itau,imisr_ztop)=
401  & fq_misr_tau_v_cth(j,itau,imisr_ztop) + boxarea
402 
403  else
404 
405  !
406  ! determine index for MISR bin set
407  !
408 
409  imisr_ztop=2
410 
411  do loop=2,n_misr_cth
412 
413  if ( box_misr_ztop(j,ibox) .gt.
414  & 1000*misr_cth_boundaries(loop+1) ) then
415 
416  imisr_ztop=loop+1
417 
418  endif
419  enddo
420 
421  if(box_cloudy(j,ibox)) then
422 
423  ! there is an isccp clouds so itau(j) is defined
424  fq_misr_tau_v_cth(j,itau,imisr_ztop)=
425  & fq_misr_tau_v_cth(j,itau,imisr_ztop) + boxarea
426 
427  else
428  ! MISR CTH resolution is trying to fill in a
429  ! broken cloud scene where there is no condensate.
430  ! The MISR CTH-1D-OD product will only put in a cloud
431  ! if the MISR cloud mask indicates cloud.
432  ! therefore we will not include this column in the histogram
433  ! in reality aerosoal and 3D effects or bright surfaces
434  ! could fool the MISR cloud mask
435 
436  ! the alternative is to count as very thin cloud ??
437 ! fq_MISR_TAU_v_CTH(1,iMISR_ztop)=
438 ! & fq_MISR_TAU_v_CTH(1,iMISR_ztop) + boxarea
439  endif
440 
441 
442  misr_mean_ztop(j)=misr_mean_ztop(j)+
443  & box_misr_ztop(j,ibox)*boxarea
444 
445  misr_cldarea(j)=misr_cldarea(j) + boxarea
446 
447  endif
448 
449  endif ! is sunlight ?
450 
451  enddo ! ibox - loop over subcolumns
452 
453  if( misr_cldarea(j) .gt. 0.) then
454  misr_mean_ztop(j)= misr_mean_ztop(j) / misr_cldarea(j) ! roj 5/2006
455  endif
456 
457  enddo ! loop over grid points
458 
459  return
460  end