source: LMDZ6/trunk/libf/phylmd/cosp2/icarus.F90 @ 3358

Last change on this file since 3358 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"

File size: 30.1 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2009, Lawrence Livemore National Security Limited Liability
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History
30! May 2015 - D. Swales - Modified for COSPv2.0
31! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32MODULE MOD_ICARUS
33  USE COSP_KINDS,          ONLY: wp
34  USE COSP_PHYS_CONSTANTS, ONLY: amd,amw,avo,grav
35  use MOD_COSP_STATS,      ONLY: hist2D
36  USE MOD_COSP_CONFIG,     ONLY: R_UNDEF,numISCCPTauBins,numISCCPPresBins,isccp_histTau, &
37                                 isccp_histPres
38  implicit none
39 
40  ! Shared Parameters                   
41  integer,parameter :: &
42       ncolprint = 0 ! Flag for debug printing (set as parameter to increase performance)
43
44  ! Cloud-top height determination
45  integer :: &
46       isccp_top_height,          & ! Top height adjustment method
47       isccp_top_height_direction   ! Direction for finding atmosphere pressure level
48
49  ! Parameters used by icarus
50  real(wp),parameter :: &
51       tauchk = -1._wp*log(0.9999999_wp), & ! Lower limit on optical depth
52       isccp_taumin = 0.3_wp,             & ! Minimum optical depth for joint-hostogram
53       pstd = 1013250._wp,                & ! Mean sea-level pressure (Pa)
54       isccp_t0 = 296._wp,                & ! Mean surface temperature (K)
55       output_missing_value = -1.E+30       ! Missing values
56
57contains
58  ! ##########################################################################
59  ! ##########################################################################
60  SUBROUTINE ICARUS(debug,debugcol,npoints,sunlit,nlev,ncol,pfull,          &
61                    phalf,qv,cc,conv,dtau_s,dtau_c,th,thd,frac_out,skt,emsfc_lw,at,&
62                    dem_s,dem_c,fq_isccp,totalcldarea, meanptop,meantaucld, &
63                    meanalbedocld, meantb,meantbclr,boxtau,boxptop,levmatch)
64       
65    ! INPUTS
66    INTEGER,intent(in) ::      & !
67         npoints,              & ! Number of model points in the horizontal
68         nlev,                 & ! Number of model levels in column
69         ncol,                 & ! Number of subcolumns
70         debug,                & ! Debug flag
71         debugcol                ! Debug column flag
72    INTEGER,intent(in),dimension(npoints) :: & !
73         sunlit                  ! 1 for day points, 0 for night time
74    REAL(WP),intent(in) ::     & !
75         emsfc_lw                ! 10.5 micron emissivity of surface (fraction) 
76    REAL(WP),intent(in),dimension(npoints) :: & !
77         skt                     ! Skin Temperature (K)
78    REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: & !
79         frac_out                ! Boxes gridbox divided up into subcolumns
80    REAL(WP),intent(in),dimension(npoints,nlev) :: & !
81         pfull,                & ! Pressure of full model levels (Pascals)
82         qv,                   & ! Water vapor specific humidity (kg vapor/ kg air)
83         cc,                   & ! Cloud cover in each model level (fraction)
84         conv,                 & ! Convective cloud cover in each model
85         at,                   & ! Temperature in each model level (K)
86         dem_c,                & ! Emissivity for convective clouds
87         dem_s,                & ! Emissivity for stratiform clouds
88         dtau_c,               & ! Optical depth for convective clouds
89         dtau_s                  ! Optical depth for stratiform clouds
90    REAL(WP),intent(in),dimension(npoints,nlev+1) :: & !
91         phalf                   ! Pressure of half model levels (Pascals)!
92    integer,intent(in) :: th,thd
93
94    ! OUTPUTS
95    REAL(WP),intent(out),dimension(npoints,7,7) :: &
96         fq_isccp                ! The fraction of the model grid box covered by clouds
97    REAL(WP),intent(out),dimension(npoints) :: &
98         totalcldarea,         & ! The fraction of model grid box columns with cloud present
99         meanptop,             & ! Mean cloud top pressure (mb) - linear averaging
100         meantaucld,           & ! Mean optical thickness
101         meanalbedocld,        & ! Mean cloud albedo 
102         meantb,               & ! Mean all-sky 10.5 micron brightness temperature
103         meantbclr               ! Mean clear-sky 10.5 micron brightness temperature
104    REAL(WP),intent(out),dimension(npoints,ncol) :: &
105         boxtau,               & ! Optical thickness in each column
106         boxptop                 ! Cloud top pressure (mb) in each column
107    INTEGER,intent(out),dimension(npoints,ncol) :: &
108         levmatch                ! Used for icarus unit testing only
109
110
111    ! INTERNAL VARIABLES
112    CHARACTER(len=10)                     :: ftn09
113    REAL(WP),dimension(npoints,ncol)      :: boxttop
114    REAL(WP),dimension(npoints,ncol,nlev) :: dtau,demIN
115    INTEGER                               :: j,ilev,ibox
116    INTEGER,dimension(nlev,ncol   )       :: acc
117
118    ! PARAMETERS
119    character ,parameter, dimension(6) :: cchar=(/' ','-','1','+','I','+'/)
120    character(len=1),parameter,dimension(6) :: cchar_realtops=(/ ' ',' ','1','1','I','I'/)
121    ! ##########################################################################
122   
123    call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dem_c,dem_s,demIN)
124    call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dtau_c,dtau_s,dtau)
125
126    call ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demIN,skt,emsfc_lw,qv,at,        &
127                      pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
128
129    call ICARUS_COLUMN(npoints,ncol,boxtau,boxptop/100._wp,sunlit,boxttop,&
130                       fq_isccp,meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
131
132    ! ##########################################################################
133    ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
134    ! ##########################################################################
135   
136    if (debugcol.ne.0) then
137       do j=1,npoints,debugcol
138         
139          ! Produce character output
140          do ilev=1,nlev
141             acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2
142             where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
143          enddo
144         
145          write(ftn09,11) j
14611        format('ftn09.',i4.4)
147          open(9, FILE=ftn09, FORM='FORMATTED')
148         
149          write(9,'(a1)') ' '
150          write(9,'(10i5)') (ilev,ilev=5,nlev,5)
151          write(9,'(a1)') ' '
152         
153          do ibox=1,ncol
154             write(9,'(40(a1),1x,40(a1))') &
155                  (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),&
156                  (cchar(acc(ilev,ibox)+1),ilev=1,nlev)
157          end do
158          close(9)
159
160       enddo       
161    end if
162   
163    return
164  end SUBROUTINE ICARUS
165 
166  ! ############################################################################
167  ! ############################################################################
168  ! ############################################################################
169  SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,at,        &
170                          pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
171    ! Inputs
172    INTEGER, intent(in) ::   &
173         ncol,               & ! Number of subcolumns
174         npoints,            & ! Number of horizontal gridpoints
175         nlev                  ! Number of vertical levels
176    INTEGER, intent(in), dimension(npoints) :: &
177         sunlit                ! 1=day 0=night
178    REAL(WP),intent(in) :: &
179         emsfc_lw              ! 10.5 micron emissivity of surface (fraction)
180    REAL(WP),intent(in), dimension(npoints) ::  &
181         skt                   ! Skin temperature
182    REAL(WP),intent(in), dimension(npoints,nlev) ::  &
183         at,                 & ! Temperature
184         pfull,              & ! Presure
185         qv                    ! Specific humidity
186    REAL(WP),intent(in), dimension(npoints,ncol,nlev) :: &
187         frac_out,           & ! Subcolumn cloud cover
188         dtau,               & ! Subcolumn optical thickness
189         demIN                 ! Subcolumn emissivity
190    REAL(WP),intent(in), dimension(npoints,nlev+1) :: &
191         phalf                 ! Pressure at model half levels
192
193    ! Outputs
194    REAL(WP),intent(inout),dimension(npoints) :: &
195         meantbclr             ! Mean clear-sky 10.5 micron brightness temperature
196    REAL(WP),intent(inout),dimension(npoints,ncol) :: &
197         boxtau,             & ! Optical thickness in each column
198         boxptop,            & ! Cloud top pressure (mb) in each column
199         boxttop               ! Cloud top temperature in each column
200    INTEGER, intent(inout),dimension(npoints,ncol)      :: levmatch
201
202    ! Local Variables
203    INTEGER :: &
204       j,ibox,ilev,k1,k2,icycle
205    INTEGER,dimension(npoints) :: &
206       nmatch,itrop
207    INTEGER,dimension(npoints,nlev-1) :: &
208       match
209    REAL(WP) :: &
210       logp,logp1,logp2,atd
211    REAL(WP),dimension(npoints) :: &
212       bb,attropmin,attrop,ptrop,atmax,btcmin,transmax,tauir,taumin,fluxtopinit,press,   &
213       dpress,atmden,rvh20,rhoave,rh20s,rfrgn,tmpexp,tauwv,wk,trans_layers_above_clrsky, &
214       fluxtop_clrsky
215    REAL(WP),dimension(npoints,nlev) :: &
216       dem_wv
217    REAL(WP),dimension(npoints,ncol) :: &
218       trans_layers_above,dem,tb,emcld,fluxtop,tau,ptop
219
220    ! ####################################################################################
221    ! Compute cloud optical depth for each column by summing up subcolumns
222    tau(1:npoints,1:ncol) = 0._wp
223    tau(1:npoints,1:ncol) = sum(dtau,dim=3)
224
225    ! Set tropopause values
226    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
227       ptrop(1:npoints)     = 5000._wp
228       attropmin(1:npoints) = 400._wp
229       atmax(1:npoints)     = 0._wp
230       attrop(1:npoints)    = 120._wp
231       itrop(1:npoints)     = 1
232
233       do ilev=1,nlev
234          where(pfull(1:npoints,ilev) .lt. 40000. .and. &
235                pfull(1:npoints,ilev) .gt.  5000. .and. &
236                at(1:npoints,ilev)    .lt. attropmin(1:npoints))
237             ptrop(1:npoints)     = pfull(1:npoints,ilev)
238             attropmin(1:npoints) = at(1:npoints,ilev)
239             attrop(1:npoints)    = attropmin(1:npoints)
240             itrop     = ilev
241          endwhere
242       enddo
243
244       do ilev=1,nlev
245          atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),&
246               at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev  .ge. itrop(1:npoints))
247       enddo
248    end if
249 
250    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
251       ! ############################################################################
252       !                        Clear-sky radiance calculation
253       !       
254       ! Compute water vapor continuum emissivity this treatment follows Schwarkzopf
255       ! and Ramasamy JGR 1999,vol 104, pages 9467-9499. The emissivity is calculated
256       ! at a wavenumber of 955 cm-1, or 10.47 microns
257       ! ############################################################################
258       do ilev=1,nlev
259          press(1:npoints)  = pfull(1:npoints,ilev)*10._wp
260          dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10
261          atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp)
262          rvh20(1:npoints)  = qv(1:npoints,ilev)*amd/amw
263          wk(1:npoints)     = rvh20(1:npoints)*avo*atmden/amd
264          rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev))
265          rh20s(1:npoints)  = rvh20(1:npoints)*rhoave(1:npoints)
266          rfrgn(1:npoints)  = rhoave(1:npoints)-rh20s(1:npoints)
267          tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0))
268          tauwv(1:npoints)  = wk(1:npoints)*1.e-20*((0.0224697_wp*rh20s(1:npoints)*      &
269                              tmpexp(1:npoints))+(3.41817e-7*rfrgn(1:npoints)))*0.98_wp
270          dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints))
271       enddo
272
273       fluxtop_clrsky(1:npoints)            = 0._wp
274       trans_layers_above_clrsky(1:npoints) = 1._wp
275       do ilev=1,nlev
276          ! Black body emission at temperature of the layer
277          bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp )
278         
279          ! Increase TOA flux by flux emitted from layer times total transmittance in layers above
280          fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
281               dem_wv(1:npoints,ilev)*bb(1:npoints)*trans_layers_above_clrsky(1:npoints)
282         
283          ! Update trans_layers_above with transmissivity from this layer for next time around loop
284          trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*&
285              (1.-dem_wv(1:npoints,ilev))               
286       enddo
287
288       ! Add in surface emission
289       bb(1:npoints) = 1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
290       fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
291           emsfc_lw * bb(1:npoints)*trans_layers_above_clrsky(1:npoints)
292
293       ! Clear Sky brightness temperature
294       meantbclr(1:npoints) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop_clrsky(1:npoints))))
295       
296       ! #################################################################################
297       !                        All-sky radiance calculation
298       ! #################################################################################
299       
300       fluxtop(1:npoints,1:ncol)            = 0._wp
301       trans_layers_above(1:npoints,1:ncol) = 1._wp
302       do ilev=1,nlev
303          ! Black body emission at temperature of the layer
304          bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp)
305         
306          do ibox=1,ncol
307             ! Emissivity
308             dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
309                                         1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
310                                         demIN(1:npoints,ibox,ilev) .eq. 0)
311
312             ! Increase TOA flux emitted from layer
313             fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)*bb*trans_layers_above(1:npoints,ibox)
314             
315             ! Update trans_layer by emitted layer from above
316             trans_layers_above(1:npoints,ibox) = trans_layers_above(1:npoints,ibox)*(1._wp-dem(1:npoints,ibox))
317          enddo
318       enddo
319
320       ! Add in surface emission
321       bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
322       do ibox=1,ncol
323          fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lw*bb(1:npoints)*trans_layers_above(1:npoints,ibox)
324       end do
325
326       ! All Sky brightness temperature
327       boxttop(1:npoints,1:ncol) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop(1:npoints,1:ncol))))
328
329       ! ################################################################################# 
330       !                            Cloud-Top Temperature
331       !
332       ! Now that you have the top of atmosphere radiance, account for ISCCP
333       ! procedures to determine cloud top temperature account for partially
334       ! transmitting cloud recompute flux ISCCP would see assuming a single layer
335       ! cloud. *NOTE* choice here of 2.13, as it is primarily ice clouds which have
336       ! partial emissivity and need the adjustment performed in this section. If it
337       ! turns out that the cloud brightness temperature is greater than 260K, then
338       ! the liquid cloud conversion factor of 2.56 is used. *NOTE* that this is
339       ! discussed on pages 85-87 of the ISCCP D level documentation
340       ! (Rossow et al. 1996)
341       ! #################################################################################
342
343       ! Compute minimum brightness temperature and optical depth
344       btcmin(1:npoints) = 1._wp /  ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp )
345
346       do ibox=1,ncol
347          transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints))
348          tauir(1:npoints)    = tau(1:npoints,ibox)/2.13_wp
349          taumin(1:npoints)   = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
350          if (isccp_top_height .eq. 1) then
351             do j=1,npoints 
352                if (transmax(j) .gt. 0.001 .and.  transmax(j) .le. 0.9999999) then
353                   fluxtopinit(j) = fluxtop(j,ibox)
354                   tauir(j) = tau(j,ibox)/2.13_wp
355                endif
356             enddo
357             do icycle=1,2
358                do j=1,npoints 
359                   if (tau(j,ibox) .gt. (tauchk)) then
360                      if (transmax(j) .gt. 0.001 .and.  transmax(j) .le. 0.9999999) then
361                         emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j)  )
362                         fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
363                         fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
364                         tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
365                         if (tb(j,ibox) .gt. 260.) then
366                            tauir(j) = tau(j,ibox) / 2.56_wp
367                         end if
368                      end if
369                   end if
370                enddo
371            enddo
372          endif
373
374          ! Cloud-top temperature
375          where(tau(1:npoints,ibox) .gt. tauchk)
376             tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox))))
377             where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints))
378                tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp
379                tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
380             endwhere
381          endwhere
382         
383          ! Clear-sky brightness temperature
384          where(tau(1:npoints,ibox) .le. tauchk)
385             tb(1:npoints,ibox) = meantbclr(1:npoints)
386          endwhere
387       enddo
388    else
389       meantbclr(1:npoints) = output_missing_value
390    end if
391
392    ! ####################################################################################
393    !                           Cloud-Top Pressure
394    !
395    ! The 2 methods differ according to whether or not you use the physical cloud
396    ! top pressure (isccp_top_height = 2) or the radiatively determined cloud top
397    ! pressure (isccp_top_height = 1 or 3)
398    ! ####################################################################################
399    do ibox=1,ncol
400       !segregate according to optical thickness
401       if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then 
402         
403          ! Find level whose temperature most closely matches brightness temperature
404          nmatch(1:npoints)=0
405          do k1=1,nlev-1
406             ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)       
407             do j=1,npoints
408                if (ilev           .ge. itrop(j)     .and. &
409                     ((at(j,ilev)  .ge. tb(j,ibox)   .and. & 
410                      at(j,ilev+1) .le. tb(j,ibox))  .or.  &
411                      (at(j,ilev)  .le. tb(j,ibox)   .and. &
412                      at(j,ilev+1) .ge. tb(j,ibox)))) then
413                   nmatch(j)=nmatch(j)+1
414                   match(j,nmatch(j))=ilev
415                endif
416             enddo
417          enddo
418
419          do j=1,npoints
420             if (nmatch(j) .ge. 1) then
421                k1 = match(j,nmatch(j))
422                k2 = k1 + 1
423                logp1 = log(pfull(j,k1))
424                logp2 = log(pfull(j,k2))
425                atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
426                logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
427                ptop(j,ibox) = exp(logp)
428                levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox)))
429             else
430                if (tb(j,ibox) .le. attrop(j)) then
431                   ptop(j,ibox)=ptrop(j)
432                   levmatch(j,ibox)=itrop(j)
433                end if
434                if (tb(j,ibox) .ge. atmax(j)) then
435                   ptop(j,ibox)=pfull(j,nlev)
436                   levmatch(j,ibox)=nlev
437                end if
438             end if
439          enddo
440       else
441          ptop(1:npoints,ibox)=0.
442          do ilev=1,nlev
443             where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0))
444                ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
445                levmatch(1:npoints,ibox)=ilev
446             endwhere
447          end do
448       end if
449       where(tau(1:npoints,ibox) .le. tauchk)
450          ptop(1:npoints,ibox)=0._wp
451          levmatch(1:npoints,ibox)=0._wp
452       endwhere
453    enddo
454
455    ! ####################################################################################
456    !                Compute subcolumn pressure and optical depth
457    ! ####################################################################################
458    boxtau(1:npoints,1:ncol)  = output_missing_value
459    boxptop(1:npoints,1:ncol) = output_missing_value
460    do ibox=1,ncol
461       do j=1,npoints
462          if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0.) then
463             if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
464                boxtau(j,ibox) = tau(j,ibox)
465                boxptop(j,ibox) = ptop(j,ibox)!/100._wp
466             endif
467          endif
468       enddo
469    enddo
470
471  end SUBROUTINE ICARUS_SUBCOLUMN
472
473  ! ######################################################################################
474  ! SUBROUTINE icarus_column
475  ! ######################################################################################
476  SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,     &
477                           meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
478    ! Inputs
479    INTEGER, intent(in) :: &
480         ncol,    & ! Number of subcolumns
481         npoints    ! Number of horizontal gridpoints
482    INTEGER, intent(in),dimension(npoints) :: &
483         sunlit     ! day=1 night=0
484    REAL(WP),intent(in),dimension(npoints,ncol) ::  &
485         boxttop,  & ! Subcolumn top temperature
486         boxptop,  & ! Subcolumn cloud top pressure
487         boxtau      ! Subcolumn optical depth
488
489    ! Outputs
490    REAL(WP),intent(inout),dimension(npoints) :: &
491         meanalbedocld, & ! Gridmean cloud albedo
492         meanptop,      & ! Gridmean cloud top pressure (mb) - linear averaging
493         meantaucld,    & ! Gridmean optical thickness
494         totalcldarea,  & ! The fraction of model grid box columns with cloud present
495         meantb           ! Gridmean all-sky 10.5 micron brightness temperature
496    REAL(WP),intent(inout),dimension(npoints,7,7) :: &
497         fq_isccp         ! The fraction of the model grid box covered by clouds
498
499    ! Local Variables
500    INTEGER :: j,ilev,ilev2
501    REAL(WP),dimension(npoints,ncol) :: albedocld
502    LOGICAL, dimension(npoints,ncol) :: box_cloudy
503
504    ! Variables for new joint-histogram implementation
505    logical,dimension(ncol) :: box_cloudy2
506
507    ! ####################################################################################
508    !                           Brightness Temperature
509    ! ####################################################################################
510    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
511       meantb(1:npoints)=sum(boxttop,2)/ncol
512    else
513       meantb(1:npoints) = output_missing_value
514    endif
515
516    ! ####################################################################################
517    !                 Determines ISCCP cloud type frequencies
518    !
519    ! Now that boxptop and boxtau have been determined, determine amount of each of the
520    ! 49 ISCCP cloud types. Also compute grid box mean cloud top pressure and
521    ! optical thickness.  The mean cloud top pressure and optical thickness are
522    ! averages over the cloudy area only. The mean cloud top pressure is a linear
523    ! average of the cloud top pressures. The mean cloud optical thickness is
524    ! computed by converting optical thickness to an albedo, averaging in albedo
525    ! units, then converting the average albedo back to a mean optical thickness. 
526    ! ####################################################################################
527
528    ! Initialize
529    albedocld(1:npoints,1:ncol)  = 0._wp
530    box_cloudy(1:npoints,1:ncol) = .false.
531   
532    ! Reset frequencies
533    !fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7)
534    do ilev=1,7
535       do ilev2=1,7
536          do j=1,npoints !
537             if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
538                fq_isccp(j,ilev,ilev2)= 0.
539             else
540                fq_isccp(j,ilev,ilev2)= output_missing_value
541             end if
542          enddo
543       enddo
544    enddo
545
546   
547    ! Reset variables need for averaging cloud properties
548    where(sunlit .eq. 1 .or. isccp_top_height .eq. 3)
549       totalcldarea(1:npoints)  = 0._wp
550       meanalbedocld(1:npoints) = 0._wp
551       meanptop(1:npoints)      = 0._wp
552       meantaucld(1:npoints)    = 0._wp
553    elsewhere
554       totalcldarea(1:npoints)  = output_missing_value
555       meanalbedocld(1:npoints) = output_missing_value
556       meanptop(1:npoints)      = output_missing_value
557       meantaucld(1:npoints)    = output_missing_value
558    endwhere
559   
560    ! Compute column quantities and joint-histogram
561    do j=1,npoints
562       ! Subcolumns that are cloudy(true) and not(false)
563       box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0.)
564
565       ! Compute joint histogram and column quantities for points that are sunlit and cloudy
566       if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then
567          ! Joint-histogram
568          call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, &
569               isccp_histPres,numISCCPPresBins,fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins))
570          fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins) = &
571               fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol
572         
573          ! Column cloud area
574          totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol
575             
576          ! Subcolumn cloud albedo
577          !albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),&
578          !     0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
579          where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
580             albedocld(j,1:ncol) = (boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp)
581          elsewhere
582             albedocld(j,1:ncol) = 0._wp
583          endwhere
584         
585          ! Column cloud albedo
586          meanalbedocld(j) = sum(albedocld(j,1:ncol))/ncol
587         
588          ! Column cloud top pressure
589          meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol
590       endif
591    enddo
592   
593    ! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
594    where(totalcldarea(1:npoints) .gt. 0)
595       meanptop(1:npoints)      = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
596       meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
597       meantaucld(1:npoints)    = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp)
598    elsewhere
599       meanptop(1:nPoints)      = output_missing_value
600       meanalbedocld(1:nPoints) = output_missing_value
601       meantaucld(1:nPoints)    = output_missing_value
602    endwhere
603    !meanptop(1:npoints)      = merge(100._wp*meanptop(1:npoints)/totalcldarea(1:npoints),&
604    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)
605    !meanalbedocld(1:npoints) = merge(meanalbedocld(1:npoints)/totalcldarea(1:npoints), &
606    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)
607    !meantaucld(1:npoints)    = merge((6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp), &
608    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)
609
610    ! Represent in percent
611    where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp
612    where(fq_isccp     .ne. output_missing_value) fq_isccp     = fq_isccp*100._wp
613   
614   
615  end SUBROUTINE ICARUS_column
616 
617  subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
618    ! INPUTS
619    integer,intent(in) :: &
620         dim1,   & ! Dimension 1 extent (Horizontal)
621         dim2,   & ! Dimension 2 extent (Subcolumn)
622         dim3      ! Dimension 3 extent (Vertical)
623    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
624         flag      ! Logical to determine the of merge var1IN and var2IN
625    real(wp),intent(in),dimension(dim1,     dim3) :: &
626         varIN1, & ! Input field 1
627         varIN2    ! Input field 2
628    ! OUTPUTS
629    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
630         varOUT    ! Merged output field
631    ! LOCAL VARIABLES
632    integer :: j
633   
634    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
635    do j=1,dim2
636       where(flag(:,j,:) .eq. 1)
637          varOUT(:,j,:) = varIN2
638       endwhere
639       where(flag(:,j,:) .eq. 2)
640          varOUT(:,j,:) = varIN1
641       endwhere
642    enddo
643  end subroutine cosp_simulator_optics
644end module MOD_ICARUS
645
Note: See TracBrowser for help on using the repository browser.