Ignore:
Timestamp:
Jul 15, 2010, 5:34:00 PM (14 years ago)
Author:
idelkadi
Message:

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90

    r1279 r1415  
    11! Copyright (c) 2009, Centre National de la Recherche Scientifique
    22! All rights reserved.
    3 ! 
    4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 
     3!
     4! Redistribution and use in source and binary forms, with or without modification, are permitted
    55! provided that the following conditions are met:
    6 ! 
    7 !     * Redistributions of source code must retain the above copyright notice, this list 
     6!
     7!     * Redistributions of source code must retain the above copyright notice, this list
    88!       of conditions and the following disclaimer.
    99!     * Redistributions in binary form must reproduce the above copyright notice, this list
    10 !       of conditions and the following disclaimer in the documentation and/or other materials 
     10!       of conditions and the following disclaimer in the documentation and/or other materials
    1111!       provided with the distribution.
    1212!     * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
    13 !       contributors may be used to endorse or promote products derived from this software without 
     13!       contributors may be used to endorse or promote products derived from this software without
    1414!       specific prior written permission.
    15 ! 
    16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
    17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
    18 ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
    19 ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
    20 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
    21 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
    22 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
     15!
     16! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
     17! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
     18! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
     19! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     20! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     21! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
     22! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
    2323! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    2424
    … …  
    3939! -----------------------------------------------------------------------------------
    4040! Lidar outputs :
    41 ! 
     41!
    4242! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction
    4343! from the lidar signals (ATB and molecular ATB) computed from model outputs
    4444!      +
    4545! Compute CFADs of lidar scattering ratio SR and of depolarization index
    46 ! 
     46!
    4747! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
    4848!
    49 ! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne : 
     49! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
    5050! - change of the cloud detection threshold S_cld from 3 to 5, for better
    5151! with both day and night observations. The optical thinest clouds are missed.
    … …  
    5353! December 2008, A. Bodas-Salcedo:
    5454! - Dimensions of pmol reduced to (npoints,llm)
     55! August 2009, A. Bodas-Salcedo:
     56! - Warning message regarding PARASOL being valid only over ocean deleted.
     57! February 2010, A. Bodas-Salcedo:
     58! - Undef passed into cosp_cfad_sr
     59! June 2010, T. Yokohata, T. Nishimura and K. Ogochi
     60! Optimisation of COSP_CFAD_SR
    5561!
    5662! Version 1.0 (June 2007)
    … …  
    7076
    7177      real undef                    ! undefined value
    72       real pnorm(npoints,ncol,llm)  ! lidar ATB 
     78      real pnorm(npoints,ncol,llm)  ! lidar ATB
    7379      real pmol(npoints,llm)        ! molecular ATB
    74       real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]   
     80      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]
    7581      real pplay(npoints,llm)       ! pressure on model levels (Pa)
    7682      logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
    … …  
    7884
    7985! c outputs :
    80       real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction 
     86      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
    8187      real cldlayer(npoints,ncat)    ! "lidar" cloud fraction (low, mid, high, total)
    82       real cfad2(npoints,max_bin,llm) ! CFADs of SR 
    83       real srbval(max_bin)           ! SR bins in CFADs 
     88      real cfad2(npoints,max_bin,llm) ! CFADs of SR
     89      real srbval(max_bin)           ! SR bins in CFADs
    8490      real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
    8591
    8692! c threshold for cloud detection :
    87       real S_clr 
    88       parameter (S_clr = 1.2) 
     93      real S_clr
     94      parameter (S_clr = 1.2)
    8995      real S_cld
    9096!      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
    … …  
    102108! c 0- Initializations
    103109! c -------------------------------------------------------
    104 ! Parasol reflectance algorithm is not valid over land. Write
    105 ! a warning if there is no land. Landmask [0 - Ocean, 1 - Land]
    106       IF ( MAXVAL(land(:)) .EQ. 0.0) THEN
    107           WRITE (*,*) 'WARNING. PARASOL reflectance is not valid over land' &
    108              & ,' and there is only land'
    109       END IF
     110!
    110111
    111112!  Should be modified in future version
    … …  
    116117! c -------------------------------------------------------
    117118!
    118 !       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
     119!       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
    119120!          x3d = pnorm/pmol
    120121!       elsewhere
    … …  
    124125      do ic = 1, ncol
    125126        pnorm_c = pnorm(:,ic,:)
    126         where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
     127        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
    127128            x3d_c = pnorm_c/pmol
    128129        elsewhere
    … …  
    142143
    143144! c -------------------------------------------------------
    144 ! c 3- CFADs 
     145! c 3- CFADs
    145146! c -------------------------------------------------------
    146147      if (ok_lidar_cfad) then
    … …  
    148149! c CFADs of subgrid-scale lidar scattering ratios :
    149150! c -------------------------------------------------------
    150       CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin, &
     151      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
    151152                 x3d, &
    152153                 S_att,S_clr,xmax,cfad2,srbval)
    … …  
    172173! if land=0 -> parasolrefl=parasolrefl
    173174        parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
    174                            + (1.0 - MAX(1.0-land(:),0.0))*undef 
     175                           + (1.0 - MAX(1.0-land(:),0.0))*undef
    175176      enddo
    176177
    177178      RETURN
    178179      END SUBROUTINE diag_lidar
    179          
    180          
     180
     181
    181182!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    182183!-------------------- FUNCTION COSP_CFAD_SR ------------------------
    183184! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
    184185!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    185       SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins, &
     186      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
    186187                      x,S_att,S_clr,xmax,cfad,srbval)
    187188      IMPLICIT NONE
    … …  
    205206! Input arguments
    206207      integer Npoints,Ncolumns,Nlevels,Nbins
    207       real xmax,S_att,S_clr
     208      real xmax,S_att,S_clr,undef
    208209! Input-output arguments
    209210      real x(Npoints,Ncolumns,Nlevels)
    … …  
    213214! Local variables
    214215      integer i, j, k, ib
     216      real srbval_ext(0:Nbins)
    215217
    216218! c -------------------------------------------------------
    … …  
    235237      cfad(:,:,:) = 0.0
    236238
    237 
     239      srbval_ext(1:Nbins) = srbval
     240      srbval_ext(0) = -1.0
    238241! c -------------------------------------------------------
    239242! c c- Compute CFAD
    240243! c -------------------------------------------------------
    241244
    242         do j = Nlevels, 1, -1
    243           do k = 1, Ncolumns
    244               where ( x(:,k,j).le.srbval(1) ) &
    245                         cfad(:,1,j) = cfad(:,1,j) + 1.0
    246           enddo  !k
    247         enddo  !j
    248 
    249       do ib = 2, Nbins
    250         do j = Nlevels, 1, -1
    251           do k = 1, Ncolumns
    252               where ( ( x(:,k,j).gt.srbval(ib-1) ) .and. ( x(:,k,j).le.srbval(ib) ) ) &
    253                         cfad(:,ib,j) = cfad(:,ib,j) + 1.0
    254           enddo  !k
    255         enddo  !j
    256       enddo  !ib
    257 
    258       cfad(:,:,:) = cfad(:,:,:) / float(Ncolumns)
     245      do j = 1, Nlevels
     246         do ib = 1, Nbins
     247            do k = 1, Ncolumns
     248               do i = 1, Npoints
     249                  if (x(i,k,j) /= undef) then
     250                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
     251                          cfad(i,ib,j) = cfad(i,ib,j) + 1.0
     252                  else
     253                     cfad(i,ib,j) = undef
     254                  endif
     255               enddo
     256            enddo
     257         enddo
     258      enddo
     259
     260      where (cfad .ne. undef)  cfad = cfad / float(Ncolumns)
    259261
    260262! c -------------------------------------------------------
    … …  
    264266!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    265267!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
    266 ! c Purpose: Cloud fraction diagnosed from lidar measurements 
     268! c Purpose: Cloud fraction diagnosed from lidar measurements
    267269!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    268270      SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, &
    … …  
    288290      real nsub(Npoints,Nlevels)
    289291
    290 
    291 ! ---------------------------------------------------------------
    292 ! 1- initialization
     292      real cldlay1(Npoints,Ncolumns)
     293      real cldlay2(Npoints,Ncolumns)
     294      real cldlay3(Npoints,Ncolumns)
     295      real nsublay1(Npoints,Ncolumns)
     296      real nsublay2(Npoints,Ncolumns)
     297      real nsublay3(Npoints,Ncolumns)
     298
     299
     300! ---------------------------------------------------------------
     301! 1- initialization
    293302! ---------------------------------------------------------------
    294303
    … …  
    317326
    318327! number of usefull sub-columns:
    319          where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  ) 
     328         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
    320329           srok(:,:,k)=1.0
    321330         elsewhere
    … …  
    329338! categories) :
    330339! ---------------------------------------------------------------
    331 ! Test abderr
     340      lidarcld = 0.0
     341      nsub = 0.0
     342
     343!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
     344      cldlay1 = 0.0
     345      cldlay2 = 0.0
     346      cldlay3 = 0.0
     347      cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
     348      nsublay1 = 0.0
     349      nsublay2 = 0.0
     350      nsublay3 = 0.0
     351      nsublay(:,:,4) = 0.0
    332352      do k = Nlevels, 1, -1
    333353       do ic = 1, Ncolumns
    334354        do ip = 1, Npoints
    335           iz=1
    336           p1 = pplay(ip,k)
    337           if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
    338             iz=3
    339           else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
    340             iz=2
     355         p1 = pplay(ip,k)
     356
     357         if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
     358            cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
     359            nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
     360         else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
     361            cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
     362            nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
     363         else
     364            cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
     365            nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
    341366         endif
    342367
    343          cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
    344          cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
     368         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
    345369         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
    346 
    347          nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
    348370         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
    349371         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
    350 
    351372        enddo
    352373       enddo
    353374      enddo
     375      cldlay(:,:,1) = cldlay1
     376      cldlay(:,:,2) = cldlay2
     377      cldlay(:,:,3) = cldlay3
     378      nsublay(:,:,1) = nsublay1
     379      nsublay(:,:,2) = nsublay2
     380      nsublay(:,:,3) = nsublay3
    354381
    355382! -- grid-box 3D cloud fraction
    … …  
    369396       do ic = 1, Ncolumns
    370397
    371           cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)   
    372           nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) 
     398          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
     399          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
    373400
    374401       enddo
    … …  
    383410      END SUBROUTINE COSP_CLDFRAC
    384411! ---------------------------------------------------------------
    385          
     412
    386413END MODULE MOD_LMD_IPSL_STATS
Note: See TracChangeset for help on using the changeset viewer.