source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/cosp/MISR_simulator.F @ 3990

Last change on this file since 3990 was 3990, checked in by millour, 8 years ago

An updated LMDZ5 (based on LMDZ rev 2786) to play with.
EM

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.