source: LMDZ4/trunk/libf/cosp/MISR_simulator.F @ 5408

Last change on this file since 5408 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

File size: 13.7 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.