Ignore:
Timestamp:
Jan 27, 2016, 10:42:32 AM (8 years ago)
Author:
idelkadi
Message:

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/phylmd/cosp/lmd_ipsl_stats.F90

    r1907 r2428  
    11! Copyright (c) 2009, Centre National de la Recherche Scientifique
    22! All rights reserved.
     3! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
     4! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lmd_ipsl_stats.F90 $
    35!
    46! Redistribution and use in source and binary forms, with or without modification, are permitted
     
    3335CONTAINS
    3436      SUBROUTINE diag_lidar(npoints,ncol,llm,max_bin,nrefl &
    35                   ,pnorm,pmol,refl,land,pplay,undef,ok_lidar_cfad &
    36                   ,cfad2,srbval &
    37                   ,ncat,lidarcld,cldlayer,parasolrefl)
     37                  ,tmp,pnorm,pnorm_perp,pmol,refl,land,pplay,undef,ok_lidar_cfad &
     38                  ,cfad2,srbval,ncat,lidarcld,lidarcldphase,cldlayer,cldlayerphase &
     39                  ,lidarcldtmp,parasolrefl)
    3840!
    3941! -----------------------------------------------------------------------------------
    4042! Lidar outputs :
    4143!
    42 ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction
    43 ! from the lidar signals (ATB and molecular ATB) computed from model outputs
     44! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction)
     45! and phase cloud fraction (3D, low/mid/high/total and 3D temperature)
     46! from the lidar signals (ATB, ATBperp and molecular ATB) computed from model outputs
    4447!      +
    4548! Compute CFADs of lidar scattering ratio SR and of depolarization index
     
    6063! Optimisation of COSP_CFAD_SR
    6164!
    62 ! Version 1.0 (June 2007)
    63 ! Version 1.1 (May 2008)
    64 ! Version 1.2 (June 2008)
    65 ! Version 2.0 (October 2008)
    66 ! Version 2.1 (December 2008)
    67 ! c------------------------------------------------------------------------------------
     65! January 2013, G. Cesana, H. Chepfer:
     66! - Add the perpendicular component of the backscattered signal (pnorm_perp) in the arguments
     67! - Add the temperature (tmp) in the arguments
     68! - Add the 3D Phase cloud fraction (lidarcldphase) in the arguments
     69! - Add the Phase low mid high cloud fraction (cldlayerphase) in the arguments
     70! - Add the 3D Phase cloud fraction as a function of temperature (lidarcldtmp) in the arguments
     71! - Modification of the phase diagnosis within the COSP_CLDFRAC routine to integrate the phase
     72!   diagnosis (3D, low/mid/high, 3D temperature)
     73! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase
     74! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376
     75!
     76! ------------------------------------------------------------------------------------
    6877
    6978! c inputs :
     
    8291      logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
    8392      real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
     93      real tmp(npoints,llm)         ! temp at each levels
     94      real pnorm_perp(npoints,ncol,llm)  ! lidar perpendicular ATB
    8495
    8596! c outputs :
    8697      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
    87       real cldlayer(npoints,ncat)    ! "lidar" cloud fraction (low, mid, high, total)
     98      real sub(npoints,llm)     ! 3D "lidar" indice
     99      real cldlayer(npoints,ncat)    ! "lidar" cloud layer fraction (low, mid, high, total)
     100
    88101      real cfad2(npoints,max_bin,llm) ! CFADs of SR
    89102      real srbval(max_bin)           ! SR bins in CFADs
     
    94107      parameter (S_clr = 1.2)
    95108      real S_cld
    96 !      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
    97       parameter (S_cld = 5.0)  ! New (dec 2008) thresold for cloud detection
     109      parameter (S_cld = 5.)  ! Thresold for cloud detection
    98110      real S_att
    99111      parameter (S_att = 0.01)
    100112
    101113! c local variables :
    102       integer ic,k
     114      integer ic,k,i,j
    103115      real x3d(npoints,ncol,llm)
    104116      real x3d_c(npoints,llm),pnorm_c(npoints,llm)
    105117      real xmax
     118
     119! Output variables
     120      integer,parameter :: nphase = 6 ! nb of cloud layer phase types (ice,liquid,undefined,false ice,false liquid,Percent of ice)
     121      real lidarcldphase(npoints,llm,nphase)   ! 3D "lidar" phase cloud fraction
     122      real lidarcldtmp(npoints,40,5)          ! 3D "lidar" phase cloud fraction as a function of temp
     123      real cldlayerphase(npoints,ncat,nphase)  ! "lidar" phase low mid high cloud fraction
     124
     125! SR detection threshold
     126      real, parameter  ::  S_cld_att = 30. ! New threshold for undefine cloud phase detection   
     127
     128
    106129!
    107130! c -------------------------------------------------------
     
    109132! c -------------------------------------------------------
    110133!
    111 
    112134!  Should be modified in future version
    113135      xmax=undef-1.0
     
    116138! c 1- Lidar scattering ratio :
    117139! c -------------------------------------------------------
    118 !
    119 !       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
    120 !          x3d = pnorm/pmol
    121 !       elsewhere
    122 !           x3d = undef
    123 !       end where
    124 ! A.B-S: pmol reduced to 2D (npoints,llm) (Dec 08)
     140
    125141      do ic = 1, ncol
    126142        pnorm_c = pnorm(:,ic,:)
     
    130146            x3d_c = undef
    131147        end where
    132         x3d(:,ic,:) = x3d_c
     148         x3d(:,ic,:) = x3d_c
    133149      enddo
    134150
     
    138154! c -------------------------------------------------------
    139155
    140       CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,  &
    141               x3d,pplay, S_att,S_cld,undef,lidarcld, &
    142               cldlayer)
     156    CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase,  &
     157              tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, &
     158              cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp)
    143159
    144160! c -------------------------------------------------------
     
    242258! c c- Compute CFAD
    243259! c -------------------------------------------------------
    244 
    245260      do j = 1, Nlevels
    246261         do ib = 1, Nbins
     
    264279      END SUBROUTINE COSP_CFAD_SR
    265280
     281
    266282!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    267283!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
    268284! c Purpose: Cloud fraction diagnosed from lidar measurements
    269285!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    270       SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, &
    271                   x,pplay,S_att,S_cld,undef,lidarcld, &
    272                   cldlayer)
     286      SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase, &
     287                  tmp,x,ATB,ATBperp,pplay,S_att,S_cld,S_cld_att,undef,lidarcld, &
     288                  cldlayer,lidarcldphase,nsub,cldlayerphase,lidarcldtemp)
     289
     290
    273291      IMPLICIT NONE
    274292! Input arguments
    275293      integer Npoints,Ncolumns,Nlevels,Ncat
    276294      real x(Npoints,Ncolumns,Nlevels)
     295
     296
     297! Local parameters
     298      integer nphase ! nb of cloud layer phase types
     299                                      ! (ice,liquid,undefined,false ice,false liquid,Percent of ice)
     300      integer,parameter  ::  Ntemp=40 ! indice of the temperature vector
     301      integer ip, k, iz, ic, ncol, nlev, i, itemp  ! loop indice
     302      real  S_cld_att ! New threshold for undefine cloud phase detection (SR=30)       
     303      integer toplvlsat  ! level of the first cloud with SR>30
     304      real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase
     305                                                                ! discrimination line   
     306
     307! Input variables
     308      real tmp(Npoints,Nlevels)                 ! temperature
     309      real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter
     310      real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter
    277311      real pplay(Npoints,Nlevels)
    278312      real S_att,S_cld
    279313      real undef
    280 ! Output :
     314
     315! Output variables
     316      real lidarcldtemp(Npoints,Ntemp,5) ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
     317      real tempmod(Ntemp+1)     ! temperature bins
     318      real lidarcldphase(Npoints,Nlevels,Nphase)    ! 3D cloud phase fraction
     319      real cldlayerphase(Npoints,Ncat,Nphase) ! low, middle, high, total cloud fractions for ice liquid and undefine phase
    281320      real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
    282321      real cldlayer(Npoints,Ncat)    ! low, middle, high, total cloud fractions
     322
    283323! Local variables
    284       integer ip, k, iz, ic
     324      real tmpi(Npoints,Ncolumns,Nlevels)       ! temperature of ice cld
     325      real tmpl(Npoints,Ncolumns,Nlevels)       ! temperature of liquid cld
     326      real tmpu(Npoints,Ncolumns,Nlevels)       ! temperature of undef cld
     327
     328      real checktemp, ATBperp_tmp ! temporary variable
     329      real checkcldlayerphase, checkcldlayerphase2 ! temporary variable
     330      real sumlidarcldtemp(Npoints,Ntemp) ! temporary variable
     331
     332      real cldlayphase(Npoints,Ncolumns,Ncat,Nphase) ! subgrided low mid high phase cloud fraction
     333      real cldlayerphasetmp(Npoints,Ncat) ! temporary variable
     334      real cldlayerphasesum(Npoints,Ncat) ! temporary variable
     335      real lidarcldtempind(Npoints,Ntemp) ! 3D Temperature indice
     336      real lidarcldphasetmp(Npoints,Nlevels)  ! 3D sum of ice and liquid cloud occurences
     337
     338
     339! Local variables
    285340      real p1
    286341      real cldy(Npoints,Ncolumns,Nlevels)
     
    290345      real nsub(Npoints,Nlevels)
    291346
     347#ifdef SYS_SX
    292348      real cldlay1(Npoints,Ncolumns)
    293349      real cldlay2(Npoints,Ncolumns)
     
    296352      real nsublay2(Npoints,Ncolumns)
    297353      real nsublay3(Npoints,Ncolumns)
     354#endif
     355
     356
    298357
    299358
     
    311370      cldlay = 0.0
    312371      nsublay = 0.0
     372
     373      ATBperp_tmp = 0.
     374      lidarcldphase(:,:,:) = 0.
     375      cldlayphase(:,:,:,:) = 0.
     376      cldlayerphase(:,:,:) = 0.
     377      tmpi(:,:,:) = 0.
     378      tmpl(:,:,:) = 0.
     379      tmpu(:,:,:) = 0.
     380      cldlayerphasesum(:,:) = 0.
     381      lidarcldtemp(:,:,:) = 0.
     382      lidarcldtempind(:,:) = 0.
     383      sumlidarcldtemp(:,:) = 0.
     384      toplvlsat=0
     385      lidarcldphasetmp(:,:) = 0.
     386
     387! temperature bins
     388      tempmod=(/-273.15,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57., &
     389                -54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18.,  &
     390                -15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /)
     391       
     392! convert C to K
     393      tempmod=tempmod+273.15
     394
     395! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
     396! (Cesana and Chepfer, JGR, 2013)
     397! ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + ATB*epsilon50 + zeta50
     398      alpha50   = 9.0322e+15
     399      beta50    = -2.1358e+12
     400      gamma50   = 173.3963e06
     401      delta50   = -3.9514e03
     402      epsilon50 = 0.2559
     403      zeta50    = -9.4776e-07
     404
    313405
    314406! ---------------------------------------------------------------
     
    334426      enddo ! k
    335427
     428
    336429! ---------------------------------------------------------------
    337430! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
     
    340433      lidarcld = 0.0
    341434      nsub = 0.0
    342 
     435#ifdef SYS_SX
    343436!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
    344437      cldlay1 = 0.0
     
    350443      nsublay3 = 0.0
    351444      nsublay(:,:,4) = 0.0
     445
    352446      do k = Nlevels, 1, -1
    353447       do ic = 1, Ncolumns
    354448        do ip = 1, Npoints
     449
     450         if(srok(ip,ic,k).gt.0.)then
     451           ! Computation of the cloud fraction as a function of the temperature
     452           ! instead of height, for ice,liquid and all clouds
     453           do itemp=1,Ntemp
     454             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
     455               lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
     456             endif
     457           enddo
     458         endif
     459
     460         if (cldy(ip,ic,k).eq.1.) then
     461           do itemp=1,Ntemp
     462             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
     463               lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
     464             endif
     465           enddo
     466         endif
     467
    355468         p1 = pplay(ip,k)
    356469
     
    379492      nsublay(:,:,2) = nsublay2
    380493      nsublay(:,:,3) = nsublay3
     494#else
     495      cldlay = 0.0
     496      nsublay = 0.0
     497      do k = Nlevels, 1, -1
     498       do ic = 1, Ncolumns
     499        do ip = 1, Npoints
     500
     501          ! Computation of the cloud fraction as a function of the temperature
     502          ! instead of height, for ice,liquid and all clouds
     503          if(srok(ip,ic,k).gt.0.)then
     504          do itemp=1,Ntemp
     505            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
     506              lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
     507            endif
     508          enddo
     509          endif
     510
     511          if(cldy(ip,ic,k).eq.1.)then
     512          do itemp=1,Ntemp
     513            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
     514              lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
     515            endif
     516          enddo
     517          endif
     518!
     519
     520          iz=1
     521          p1 = pplay(ip,k)
     522          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
     523            iz=3
     524          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
     525            iz=2
     526         endif
     527
     528         cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
     529         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
     530         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
     531
     532         nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
     533         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
     534         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
     535
     536        enddo
     537       enddo
     538      enddo
     539#endif
     540
    381541
    382542! -- grid-box 3D cloud fraction
     
    407567      endwhere
    408568
    409       RETURN
     569! ---------------------------------------------------------------
     570! 4- grid-box 3D cloud Phase :
     571! ---------------------------------------------------------------
     572! ---------------------------------------------------------------
     573! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km
     574! ---------------------------------------------------------------
     575do ncol=1,Ncolumns
     576do i=1,Npoints
     577
     578      do nlev=Nlevels,18,-1  ! from 19.2km until 8.16km
     579         p1 = pplay(i,nlev)
     580
     581
     582! Avoid zero values
     583        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
     584! Computation of the ATBperp along the phase discrimination line
     585           ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
     586                         (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
     587                          ATB(i,ncol,nlev)*epsilon50 + zeta50
     588
     589!____________________________________________________________________________________________________
     590!
     591!4.1.a Ice: ATBperp above the phase discrimination line
     592!____________________________________________________________________________________________________
     593!
     594           if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
     595             ! ICE with temperature above 273,15°K = Liquid (false ice)
     596            if(tmp(i,nlev).gt.273.15)then                ! Temperature above 273,15 K
     597              ! Liquid: False ice corrected by the temperature to Liquid
     598               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.   ! false ice detection ==> added to Liquid
     599               tmpl(i,ncol,nlev)=tmp(i,nlev)
     600               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.   ! keep the information "temperature criterium used"
     601                                                    ! to classify the phase cloud
     602                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     603                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     604                   cldlayphase(i,ncol,3,2) = 1.
     605                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     606                   cldlayphase(i,ncol,2,2) = 1.
     607                else                                                    ! low cloud
     608                   cldlayphase(i,ncol,1,2) = 1.
     609                endif
     610                   cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
     611                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     612                   cldlayphase(i,ncol,3,5) = 1.
     613                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     614                   cldlayphase(i,ncol,2,5) = 1.
     615                else                                                    ! low cloud
     616                   cldlayphase(i,ncol,1,5) = 1.
     617                endif
     618
     619             else
     620             ! ICE with temperature below 273,15°K
     621              lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
     622              tmpi(i,ncol,nlev)=tmp(i,nlev)
     623                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     624                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     625                   cldlayphase(i,ncol,3,1) = 1.
     626                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     627                   cldlayphase(i,ncol,2,1) = 1.
     628                else                                                    ! low cloud
     629                   cldlayphase(i,ncol,1,1) = 1.
     630                endif
     631
     632              endif
     633
     634!____________________________________________________________________________________________________
     635!
     636! 4.1.b Liquid: ATBperp below the phase discrimination line
     637!____________________________________________________________________________________________________
     638!
     639             else                                        ! Liquid clouds
     640              ! Liquid with temperature above 231,15°K
     641            if(tmp(i,nlev).gt.231.15)then
     642               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
     643               tmpl(i,ncol,nlev)=tmp(i,nlev)
     644                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     645                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     646                   cldlayphase(i,ncol,3,2) = 1. 
     647                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     648                   cldlayphase(i,ncol,2,2) = 1.
     649                else                                                    ! low cloud
     650                   cldlayphase(i,ncol,1,2) = 1.
     651                endif
     652
     653             else
     654             ! Liquid with temperature below 231,15°K = Ice (false liquid)
     655               tmpi(i,ncol,nlev)=tmp(i,nlev)
     656               lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.   ! false liquid detection ==> added to ice
     657               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.   ! keep the information "temperature criterium used"
     658                                                    ! to classify the phase cloud
     659                   cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
     660                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     661                   cldlayphase(i,ncol,3,4) = 1. 
     662                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     663                   cldlayphase(i,ncol,2,4) = 1.
     664                else                                                    ! low cloud
     665                   cldlayphase(i,ncol,1,4) = 1.
     666                endif
     667                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     668                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     669                   cldlayphase(i,ncol,3,1) = 1. 
     670                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     671                   cldlayphase(i,ncol,2,1) = 1.
     672                else                                                    ! low cloud
     673                   cldlayphase(i,ncol,1,1) = 1.
     674                endif
     675
     676             endif
     677
     678            endif  ! end of discrimination condition
     679         endif  ! end of cloud condition
     680      enddo ! end of altitude loop
     681
     682
     683
     684! ---------------------------------------------------------------
     685! 4.2 - For Cloudy pixels with 0km < z < 8.16km
     686! ---------------------------------------------------------------
     687
     688      toplvlsat=0
     689      do nlev=17,1,-1  ! from 8.16km until 0km
     690         p1 = pplay(i,nlev)
     691
     692        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
     693! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50
     694!                                  + ATB*epsilon50 + zeta50
     695! Computation of the ATBperp of the phase discrimination line
     696           ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
     697                         (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
     698                          ATB(i,ncol,nlev)*epsilon50 + zeta50
     699!____________________________________________________________________________________________________
     700!
     701! 4.2.a Ice: ATBperp above the phase discrimination line
     702!____________________________________________________________________________________________________
     703!
     704            ! ICE with temperature above 273,15°K = Liquid (false ice)
     705          if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
     706            if(tmp(i,nlev).gt.273.15)then
     707               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.  ! false ice ==> liq
     708               tmpl(i,ncol,nlev)=tmp(i,nlev)
     709               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.
     710
     711                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     712               if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
     713                   cldlayphase(i,ncol,3,2) = 1.
     714                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     715                   cldlayphase(i,ncol,2,2) = 1.
     716                else                                                    ! low cloud
     717                   cldlayphase(i,ncol,1,2) = 1.
     718                endif
     719
     720                   cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
     721                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     722                   cldlayphase(i,ncol,3,5) = 1.
     723                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     724                   cldlayphase(i,ncol,2,5) = 1.
     725                else                                                    ! low cloud
     726                   cldlayphase(i,ncol,1,5) = 1.
     727                endif
     728
     729             else
     730              ! ICE with temperature below 273,15°K
     731              lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
     732              tmpi(i,ncol,nlev)=tmp(i,nlev)
     733
     734                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     735                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     736                   cldlayphase(i,ncol,3,1) = 1.
     737                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     738                   cldlayphase(i,ncol,2,1) = 1.
     739                else                                                    ! low cloud
     740                   cldlayphase(i,ncol,1,1) = 1.
     741                endif
     742
     743              endif
     744
     745!____________________________________________________________________________________________________
     746!
     747! 4.2.b Liquid: ATBperp below the phase discrimination line
     748!____________________________________________________________________________________________________
     749!
     750          else 
     751             ! Liquid with temperature above 231,15°K
     752            if(tmp(i,nlev).gt.231.15)then
     753               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
     754               tmpl(i,ncol,nlev)=tmp(i,nlev)
     755
     756                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     757                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     758                   cldlayphase(i,ncol,3,2) = 1. 
     759                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     760                   cldlayphase(i,ncol,2,2) = 1.
     761                else                                                    ! low cloud
     762                   cldlayphase(i,ncol,1,2) = 1.
     763                endif
     764
     765             else
     766             ! Liquid with temperature below 231,15°K = Ice (false liquid)
     767               tmpi(i,ncol,nlev)=tmp(i,nlev)
     768               lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.  ! false liq ==> ice
     769               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.  ! false liq ==> ice
     770
     771                   cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
     772                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     773                   cldlayphase(i,ncol,3,4) = 1. 
     774                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     775                   cldlayphase(i,ncol,2,4) = 1.
     776                else                                                    ! low cloud
     777                   cldlayphase(i,ncol,1,4) = 1.
     778                endif
     779
     780                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     781                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     782                   cldlayphase(i,ncol,3,1) = 1. 
     783                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     784                   cldlayphase(i,ncol,2,1) = 1.
     785                else                                                    ! low cloud
     786                   cldlayphase(i,ncol,1,1) = 1.
     787                endif
     788
     789             endif
     790           endif  ! end of discrimination condition
     791
     792            toplvlsat=0
     793
     794           ! Find the level of the highest cloud with SR>30
     795            if(x(i,ncol,nlev).gt.S_cld_att)then  ! SR > 30.
     796                toplvlsat=nlev-1
     797                goto 99
     798            endif
     799
     800        endif  ! end of cloud condition
     801       enddo  ! end of altitude loop
     802
     80399 continue
     804
     805!____________________________________________________________________________________________________
     806!
     807! Undefined phase: For a cloud located below another cloud with SR>30
     808! see Cesana and Chepfer 2013 Sect.III.2
     809!____________________________________________________________________________________________________
     810!
     811if(toplvlsat.ne.0)then         
     812      do nlev=toplvlsat,1,-1
     813         p1 = pplay(i,nlev)
     814        if(cldy(i,ncol,nlev).eq.1.)then
     815           lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
     816           tmpu(i,ncol,nlev)=tmp(i,nlev)
     817
     818                   cldlayphase(i,ncol,4,3) = 1.                         ! tot cloud
     819          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
     820             cldlayphase(i,ncol,3,3) = 1.
     821          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid cloud
     822             cldlayphase(i,ncol,2,3) = 1.
     823          else                                                     ! low cloud
     824             cldlayphase(i,ncol,1,3) = 1.
     825          endif
     826
     827        endif   
     828      enddo
     829endif
     830     
     831      toplvlsat=0
     832
     833enddo
     834enddo
     835
     836
     837
     838!____________________________________________________________________________________________________
     839!
     840! Computation of final cloud phase diagnosis
     841!____________________________________________________________________________________________________
     842!
     843
     844! Compute the Ice percentage in cloud = ice/(ice+liq) as a function
     845! of the occurrences
     846lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
     847WHERE (lidarcldphasetmp(:,:).gt. 0.)
     848   lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
     849ELSEWHERE
     850   lidarcldphase(:,:,6) = undef
     851ENDWHERE
     852
     853! Compute Phase 3D Cloud Fraction
     854     WHERE ( nsub(:,:).gt.0.0 )
     855       lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
     856       lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
     857       lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
     858       lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
     859       lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
     860     ELSEWHERE
     861       lidarcldphase(:,:,1) = undef
     862       lidarcldphase(:,:,2) = undef
     863       lidarcldphase(:,:,3) = undef
     864       lidarcldphase(:,:,4) = undef
     865       lidarcldphase(:,:,5) = undef
     866     ENDWHERE
     867
     868
     869! Compute Phase low mid high cloud fractions
     870    do iz = 1, Ncat
     871       do i=1,Nphase-3
     872       do ic = 1, Ncolumns
     873          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
     874          cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i)
     875       enddo
     876      enddo
     877    enddo
     878
     879    do iz = 1, Ncat
     880       do i=4,5
     881       do ic = 1, Ncolumns
     882          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)         
     883       enddo
     884       enddo
     885    enddo
     886   
     887! Compute the Ice percentage in cloud = ice/(ice+liq)
     888cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
     889    WHERE (cldlayerphasetmp(:,:).gt. 0.)
     890       cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
     891    ELSEWHERE
     892       cldlayerphase(:,:,6) = undef
     893    ENDWHERE
     894
     895    do i=1,Nphase-1
     896      WHERE ( cldlayerphasesum(:,:).gt.0.0 )
     897         cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
     898      ENDWHERE
     899    enddo
     900
     901
     902    do i=1,Npoints
     903       do iz=1,Ncat
     904          checkcldlayerphase=0.
     905          checkcldlayerphase2=0.
     906
     907          if (cldlayerphasesum(i,iz).gt.0.0 )then
     908             do ic=1,Nphase-3
     909                checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic) 
     910             enddo
     911             checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase
     912             if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
     913
     914          endif
     915
     916       enddo
     917    enddo
     918
     919    do i=1,Nphase-1
     920      WHERE ( nsublayer(:,:).eq.0.0 )
     921         cldlayerphase(:,:,i) = undef
     922      ENDWHERE
     923   enddo
     924
     925
     926
     927! Compute Phase 3D as a function of temperature
     928do nlev=1,Nlevels
     929do ncol=1,Ncolumns     
     930do i=1,Npoints
     931do itemp=1,Ntemp
     932if(tmpi(i,ncol,nlev).gt.0.)then
     933      if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
     934        lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1.
     935      endif
     936elseif(tmpl(i,ncol,nlev).gt.0.)then
     937      if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then
     938        lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1.
     939      endif
     940elseif(tmpu(i,ncol,nlev).gt.0.)then
     941      if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then
     942        lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1.
     943      endif
     944endif
     945enddo
     946enddo
     947enddo
     948enddo
     949
     950! Check temperature cloud fraction
     951do i=1,Npoints
     952   do itemp=1,Ntemp
     953checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
     954
     955        if(checktemp.NE.lidarcldtemp(i,itemp,1))then
     956          print *, i,itemp
     957          print *, lidarcldtemp(i,itemp,1:4)
     958        endif
     959
     960   enddo
     961enddo
     962
     963! Compute the Ice percentage in cloud = ice/(ice+liq)
     964!   sumlidarcldtemp=sum(lidarcldtemp(:,:,2:3),3)
     965   sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)
     966
     967WHERE(sumlidarcldtemp(:,:)>0.)
     968  lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
     969ELSEWHERE
     970  lidarcldtemp(:,:,5)=undef
     971ENDWHERE
     972
     973do i=1,4
     974  WHERE(lidarcldtempind(:,:).gt.0.)
     975     lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
     976  ELSEWHERE
     977     lidarcldtemp(:,:,i) = undef
     978  ENDWHERE
     979enddo
     980
     981       RETURN
    410982      END SUBROUTINE COSP_CLDFRAC
    411983! ---------------------------------------------------------------
    412984
     985
    413986END MODULE MOD_LMD_IPSL_STATS
Note: See TracChangeset for help on using the changeset viewer.