source: LMDZ6/branches/contrails/libf/phylmd/cosp/MISR_simulator.f90 @ 5456

Last change on this file since 5456 was 5248, checked in by abarral, 2 months ago

(cont.) Convert fixed-form to free-form sources .F -> .{f,F}90

  • 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: 14.3 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
26SUBROUTINE 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        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
475end subroutine misr_simulator
Note: See TracBrowser for help on using the repository browser.