! Copyright (c) 2009, Centre National de la Recherche Scientifique ! All rights reserved. ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lmd_ipsl_stats.F90 $ ! ! Redistribution and use in source and binary forms, with or without modification, are permitted ! provided that the following conditions are met: ! ! * Redistributions of source code must retain the above copyright notice, this list ! of conditions and the following disclaimer. ! * Redistributions in binary form must reproduce the above copyright notice, this list ! of conditions and the following disclaimer in the documentation and/or other materials ! provided with the distribution. ! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its ! contributors may be used to endorse or promote products derived from this software without ! specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. !------------------------------------------------------------------------------------ ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France). !------------------------------------------------------------------------------------ MODULE MOD_LMD_IPSL_STATS USE MOD_LLNL_STATS IMPLICIT NONE CONTAINS SUBROUTINE diag_lidar(npoints,ncol,llm,max_bin,nrefl & ,tmp,pnorm,pnorm_perp,pmol,refl,land,pplay,undef,ok_lidar_cfad & ,cfad2,srbval,ncat,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ ,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR) !OPAQ !TIBO ! ! ----------------------------------------------------------------------------------- ! Lidar outputs : ! ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction) ! and phase cloud fraction (3D, low/mid/high/total and 3D temperature) ! from the lidar signals (ATB, ATBperp and molecular ATB) computed from model outputs ! + ! Compute CFADs of lidar scattering ratio SR and of depolarization index ! ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France). ! ! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : ! - change of the cloud detection threshold S_cld from 3 to 5, for better ! with both day and night observations. The optical thinest clouds are missed. ! - remove of the detection of the first fully attenuated layer encountered from above. ! December 2008, A. Bodas-Salcedo: ! - Dimensions of pmol reduced to (npoints,llm) ! August 2009, A. Bodas-Salcedo: ! - Warning message regarding PARASOL being valid only over ocean deleted. ! February 2010, A. Bodas-Salcedo: ! - Undef passed into cosp_cfad_sr ! June 2010, T. Yokohata, T. Nishimura and K. Ogochi ! Optimisation of COSP_CFAD_SR ! ! January 2013, G. Cesana, H. Chepfer: ! - Add the perpendicular component of the backscattered signal (pnorm_perp) in the arguments ! - Add the temperature (tmp) in the arguments ! - Add the 3D Phase cloud fraction (lidarcldphase) in the arguments ! - Add the Phase low mid high cloud fraction (cldlayerphase) in the arguments ! - Add the 3D Phase cloud fraction as a function of temperature (lidarcldtmp) in the arguments ! - Modification of the phase diagnosis within the COSP_CLDFRAC routine to integrate the phase ! diagnosis (3D, low/mid/high, 3D temperature) ! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase ! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376 ! ! ------------------------------------------------------------------------------------ ! c inputs : integer npoints integer ncol integer llm integer max_bin ! nb of bins for SR CFADs integer ncat ! nb of cloud layer types (low,mid,high,total) integer ntype ! nb of OPAQ products (opaque and thin clouds, z_opaque) !OPAQ integer nrefl ! nb of solar zenith angles for parasol reflectances real undef ! undefined value real pnorm(npoints,ncol,llm) ! lidar ATB real pmol(npoints,llm) ! molecular ATB real land(npoints) ! Landmask [0 - Ocean, 1 - Land] real pplay(npoints,llm) ! pressure on model levels (Pa) logical ok_lidar_cfad ! true if lidar CFAD diagnostics need to be computed real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol real tmp(npoints,llm) ! temp at each levels real pnorm_perp(npoints,ncol,llm) ! lidar perpendicular ATB real vgrid_z(llm) ! mid-level altitude of the output vertical grid !OPAQ ! c outputs : real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction real lidarcldtype(npoints,llm,ntype+1) ! 3D "lidar" OPAQ type fraction + opacity !OPAQ real sub(npoints,llm) ! 3D "lidar" indice real cldlayer(npoints,ncat) ! "lidar" cloud layer fraction (low, mid, high, total) real cldtype(npoints,ntype) ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque) !OPAQ real cfad2(npoints,max_bin,llm) ! CFADs of SR real srbval(max_bin) ! SR bins in CFADs real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance ! real profSR(npoints,ncol,llm) ! tableau avec les subcolumns SR !TIBO real profSR(npoints,llm,ncol) ! tableau avec les subcolumns SR !TIBO2 ! c threshold for cloud detection : real S_clr parameter (S_clr = 1.2) real S_cld parameter (S_cld = 5.) ! Thresold for cloud detection real S_att parameter (S_att = 0.01) ! parameter (S_att = 0.06) !OPAQ ! Threshold for "surface detection" equivalent ! c local variables : integer ic,k,i,j real x3d(npoints,ncol,llm) real x3d_c(npoints,llm),pnorm_c(npoints,llm) real xmax ! Output variables integer,parameter :: nphase = 6 ! nb of cloud layer phase types (ice,liquid,undefined,false ice,false liquid,Percent of ice) real lidarcldphase(npoints,llm,nphase) ! 3D "lidar" phase cloud fraction real lidarcldtmp(npoints,40,5) ! 3D "lidar" phase cloud fraction as a function of temp real cldlayerphase(npoints,ncat,nphase) ! "lidar" phase low mid high cloud fraction ! SR detection threshold real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection ! ! c ------------------------------------------------------- ! c 0- Initializations ! c ------------------------------------------------------- ! ! Should be modified in future version xmax=undef-1.0 ! c ------------------------------------------------------- ! c 1- Lidar scattering ratio : ! c ------------------------------------------------------- do ic = 1, ncol pnorm_c = pnorm(:,ic,:) where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) x3d_c = pnorm_c/pmol elsewhere x3d_c = undef end where x3d(:,ic,:) = x3d_c ! profSR(:,ic,:) = x3d(:,ic,:) !TIBO profSR(:,:,ic) = x3d(:,ic,:) !TIBO2 enddo ! c ------------------------------------------------------- ! c 2- Diagnose cloud fractions (3D, low, middle, high, total) ! c from subgrid-scale lidar scattering ratios : ! c ------------------------------------------------------- CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase, & tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, & cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp) CALL COSP_OPAQ(npoints,ncol,llm,ntype,x3d,S_cld,undef,lidarcldtype, & !OPAQ cldtype,vgrid_z) !OPAQ ! c ------------------------------------------------------- ! c 3- CFADs ! c ------------------------------------------------------- if (ok_lidar_cfad) then ! ! c CFADs of subgrid-scale lidar scattering ratios : ! c ------------------------------------------------------- CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, & x3d, & S_att,S_clr,xmax,cfad2,srbval) endif ! ok_lidar_cfad ! c ------------------------------------------------------- ! c ------------------------------------------------------- ! c 4- Compute grid-box averaged Parasol reflectances ! c ------------------------------------------------------- parasolrefl(:,:) = 0.0 do k = 1, nrefl do ic = 1, ncol parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k) enddo enddo do k = 1, nrefl parasolrefl(:,k) = parasolrefl(:,k) / float(ncol) ! if land=1 -> parasolrefl=undef ! if land=0 -> parasolrefl=parasolrefl parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) & + (1.0 - MAX(1.0-land(:),0.0))*undef enddo RETURN END SUBROUTINE diag_lidar !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------- FUNCTION COSP_CFAD_SR ------------------------ ! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, & x,S_att,S_clr,xmax,cfad,srbval) IMPLICIT NONE !--- Input arguments ! Npoints: Number of horizontal points ! Ncolumns: Number of subcolumns ! Nlevels: Number of levels ! Nbins: Number of x axis bins ! xmax: maximum value allowed for x ! S_att: Threshold for full attenuation ! S_clr: Threshold for clear-sky layer ! !--- Input-Outout arguments ! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs ! ! -- Output arguments ! srbval : values of the histogram bins ! cfad: 2D histogram on each horizontal point ! Input arguments integer Npoints,Ncolumns,Nlevels,Nbins real xmax,S_att,S_clr,undef ! Input-output arguments real x(Npoints,Ncolumns,Nlevels) ! Output : real cfad(Npoints,Nbins,Nlevels) real srbval(Nbins) ! Local variables integer i, j, k, ib real srbval_ext(0:Nbins) ! c ------------------------------------------------------- ! c 0- Initializations ! c ------------------------------------------------------- if ( Nbins .lt. 6) return srbval(1) = S_att srbval(2) = S_clr srbval(3) = 3.0 srbval(4) = 5.0 srbval(5) = 7.0 srbval(6) = 10.0 do i = 7, MIN(10,Nbins) srbval(i) = srbval(i-1) + 5.0 enddo DO i = 11, MIN(13,Nbins) srbval(i) = srbval(i-1) + 10.0 enddo srbval(MIN(14,Nbins)) = 80.0 srbval(Nbins) = xmax cfad(:,:,:) = 0.0 srbval_ext(1:Nbins) = srbval srbval_ext(0) = -1.0 ! c ------------------------------------------------------- ! c c- Compute CFAD ! c ------------------------------------------------------- do j = 1, Nlevels do ib = 1, Nbins do k = 1, Ncolumns do i = 1, Npoints if (x(i,k,j) /= undef) then if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) & cfad(i,ib,j) = cfad(i,ib,j) + 1.0 else cfad(i,ib,j) = undef endif enddo enddo enddo enddo where (cfad .ne. undef) cfad = cfad / float(Ncolumns) ! c ------------------------------------------------------- RETURN END SUBROUTINE COSP_CFAD_SR !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------- SUBROUTINE COSP_CLDFRAC ------------------- ! c Purpose: Cloud fraction diagnosed from lidar measurements !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase, & tmp,x,ATB,ATBperp,pplay,S_att,S_cld,S_cld_att,undef,lidarcld, & cldlayer,lidarcldphase,nsub,cldlayerphase,lidarcldtemp) IMPLICIT NONE ! Input arguments integer Npoints,Ncolumns,Nlevels,Ncat real x(Npoints,Ncolumns,Nlevels) ! Local parameters integer nphase ! nb of cloud layer phase types ! (ice,liquid,undefined,false ice,false liquid,Percent of ice) integer,parameter :: Ntemp=40 ! indice of the temperature vector integer ip, k, iz, ic, ncol, nlev, i, itemp ! loop indice real S_cld_att ! New threshold for undefine cloud phase detection (SR=30) integer toplvlsat ! level of the first cloud with SR>30 real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase ! discrimination line ! Input variables real tmp(Npoints,Nlevels) ! temperature real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter real pplay(Npoints,Nlevels) real S_att,S_cld real undef ! Output variables real lidarcldtemp(Npoints,Ntemp,5) ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq real tempmod(Ntemp+1) ! temperature bins real lidarcldphase(Npoints,Nlevels,Nphase) ! 3D cloud phase fraction real cldlayerphase(Npoints,Ncat,Nphase) ! low, middle, high, total cloud fractions for ice liquid and undefine phase real lidarcld(Npoints,Nlevels) ! 3D cloud fraction real cldlayer(Npoints,Ncat) ! low, middle, high, total cloud fractions ! Local variables real tmpi(Npoints,Ncolumns,Nlevels) ! temperature of ice cld real tmpl(Npoints,Ncolumns,Nlevels) ! temperature of liquid cld real tmpu(Npoints,Ncolumns,Nlevels) ! temperature of undef cld real checktemp, ATBperp_tmp ! temporary variable real checkcldlayerphase, checkcldlayerphase2 ! temporary variable real sumlidarcldtemp(Npoints,Ntemp) ! temporary variable real cldlayphase(Npoints,Ncolumns,Ncat,Nphase) ! subgrided low mid high phase cloud fraction real cldlayerphasetmp(Npoints,Ncat) ! temporary variable real cldlayerphasesum(Npoints,Ncat) ! temporary variable real lidarcldtempind(Npoints,Ntemp) ! 3D Temperature indice real lidarcldphasetmp(Npoints,Nlevels) ! 3D sum of ice and liquid cloud occurences ! Local variables real p1 real cldy(Npoints,Ncolumns,Nlevels) real srok(Npoints,Ncolumns,Nlevels) real cldlay(Npoints,Ncolumns,Ncat) real nsublay(Npoints,Ncolumns,Ncat), nsublayer(Npoints,Ncat) real nsub(Npoints,Nlevels) #ifdef SYS_SX real cldlay1(Npoints,Ncolumns) real cldlay2(Npoints,Ncolumns) real cldlay3(Npoints,Ncolumns) real nsublay1(Npoints,Ncolumns) real nsublay2(Npoints,Ncolumns) real nsublay3(Npoints,Ncolumns) #endif ! --------------------------------------------------------------- ! 1- initialization ! --------------------------------------------------------------- if ( Ncat .ne. 4 ) then print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat stop endif lidarcld = 0.0 nsub = 0.0 cldlay = 0.0 nsublay = 0.0 ATBperp_tmp = 0. lidarcldphase(:,:,:) = 0. cldlayphase(:,:,:,:) = 0. cldlayerphase(:,:,:) = 0. tmpi(:,:,:) = 0. tmpl(:,:,:) = 0. tmpu(:,:,:) = 0. cldlayerphasesum(:,:) = 0. lidarcldtemp(:,:,:) = 0. lidarcldtempind(:,:) = 0. sumlidarcldtemp(:,:) = 0. toplvlsat=0 lidarcldphasetmp(:,:) = 0. ! temperature bins tempmod=(/-273.15,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57., & -54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18., & -15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /) ! convert C to K tempmod=tempmod+273.15 ! Polynomial coefficient of the phase discrimination line used to separate liquid from ice ! (Cesana and Chepfer, JGR, 2013) ! ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + ATB*epsilon50 + zeta50 alpha50 = 9.0322e+15 beta50 = -2.1358e+12 gamma50 = 173.3963e06 delta50 = -3.9514e03 epsilon50 = 0.2559 zeta50 = -9.4776e-07 ! --------------------------------------------------------------- ! 2- Cloud detection ! --------------------------------------------------------------- do k = 1, Nlevels ! cloud detection at subgrid-scale: where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) ) cldy(:,:,k)=1.0 elsewhere cldy(:,:,k)=0.0 endwhere ! number of usefull sub-columns: where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) ) srok(:,:,k)=1.0 elsewhere srok(:,:,k)=0.0 endwhere enddo ! k ! --------------------------------------------------------------- ! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure ! categories) : ! --------------------------------------------------------------- lidarcld = 0.0 nsub = 0.0 #ifdef SYS_SX !! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts. cldlay1 = 0.0 cldlay2 = 0.0 cldlay3 = 0.0 cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4 nsublay1 = 0.0 nsublay2 = 0.0 nsublay3 = 0.0 nsublay(:,:,4) = 0.0 do k = Nlevels, 1, -1 do ic = 1, Ncolumns do ip = 1, Npoints if(srok(ip,ic,k).gt.0.)then ! Computation of the cloud fraction as a function of the temperature ! instead of height, for ice,liquid and all clouds do itemp=1,Ntemp if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1. endif enddo endif if (cldy(ip,ic,k).eq.1.) then do itemp=1,Ntemp if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1. endif enddo endif p1 = pplay(ip,k) if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k)) nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k)) else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid clouds cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k)) nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k)) else cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k)) nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k)) endif cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k)) lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k) nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k)) nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k) enddo enddo enddo cldlay(:,:,1) = cldlay1 cldlay(:,:,2) = cldlay2 cldlay(:,:,3) = cldlay3 nsublay(:,:,1) = nsublay1 nsublay(:,:,2) = nsublay2 nsublay(:,:,3) = nsublay3 #else cldlay = 0.0 nsublay = 0.0 do k = Nlevels, 1, -1 do ic = 1, Ncolumns do ip = 1, Npoints ! Computation of the cloud fraction as a function of the temperature ! instead of height, for ice,liquid and all clouds if(srok(ip,ic,k).gt.0.)then do itemp=1,Ntemp if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1. endif enddo endif if(cldy(ip,ic,k).eq.1.)then do itemp=1,Ntemp if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1. endif enddo endif ! iz=1 p1 = pplay(ip,k) if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds iz=3 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid clouds iz=2 endif cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k)) cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k)) lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k) nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k)) nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k)) nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k) enddo enddo enddo #endif ! -- grid-box 3D cloud fraction where ( nsub(:,:).gt.0.0 ) lidarcld(:,:) = lidarcld(:,:)/nsub(:,:) elsewhere lidarcld(:,:) = undef endwhere ! -- layered cloud fractions cldlayer = 0.0 nsublayer = 0.0 do iz = 1, Ncat do ic = 1, Ncolumns cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz) nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) enddo enddo where ( nsublayer(:,:).gt.0.0 ) cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:) elsewhere cldlayer(:,:) = undef endwhere ! --------------------------------------------------------------- ! 4- grid-box 3D cloud Phase : ! --------------------------------------------------------------- ! --------------------------------------------------------------- ! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km ! --------------------------------------------------------------- do ncol=1,Ncolumns do i=1,Npoints do nlev=Nlevels,18,-1 ! from 19.2km until 8.16km p1 = pplay(i,nlev) ! Avoid zero values if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then ! Computation of the ATBperp along the phase discrimination line ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + & (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + & ATB(i,ncol,nlev)*epsilon50 + zeta50 !____________________________________________________________________________________________________ ! !4.1.a Ice: ATBperp above the phase discrimination line !____________________________________________________________________________________________________ ! if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds ! ICE with temperature above 273,15°K = Liquid (false ice) if(tmp(i,nlev).gt.273.15)then ! Temperature above 273,15 K ! Liquid: False ice corrected by the temperature to Liquid lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice detection ==> added to Liquid tmpl(i,ncol,nlev)=tmp(i,nlev) lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. ! keep the information "temperature criterium used" ! to classify the phase cloud cldlayphase(i,ncol,4,2) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,2) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,2) = 1. else ! low cloud cldlayphase(i,ncol,1,2) = 1. endif cldlayphase(i,ncol,4,5) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,5) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,5) = 1. else ! low cloud cldlayphase(i,ncol,1,5) = 1. endif else ! ICE with temperature below 273,15°K lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. tmpi(i,ncol,nlev)=tmp(i,nlev) cldlayphase(i,ncol,4,1) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,1) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,1) = 1. else ! low cloud cldlayphase(i,ncol,1,1) = 1. endif endif !____________________________________________________________________________________________________ ! ! 4.1.b Liquid: ATBperp below the phase discrimination line !____________________________________________________________________________________________________ ! else ! Liquid clouds ! Liquid with temperature above 231,15°K if(tmp(i,nlev).gt.231.15)then lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. tmpl(i,ncol,nlev)=tmp(i,nlev) cldlayphase(i,ncol,4,2) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,2) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,2) = 1. else ! low cloud cldlayphase(i,ncol,1,2) = 1. endif else ! Liquid with temperature below 231,15°K = Ice (false liquid) tmpi(i,ncol,nlev)=tmp(i,nlev) lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liquid detection ==> added to ice lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! keep the information "temperature criterium used" ! to classify the phase cloud cldlayphase(i,ncol,4,4) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,4) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,4) = 1. else ! low cloud cldlayphase(i,ncol,1,4) = 1. endif cldlayphase(i,ncol,4,1) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,1) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,1) = 1. else ! low cloud cldlayphase(i,ncol,1,1) = 1. endif endif endif ! end of discrimination condition endif ! end of cloud condition enddo ! end of altitude loop ! --------------------------------------------------------------- ! 4.2 - For Cloudy pixels with 0km < z < 8.16km ! --------------------------------------------------------------- toplvlsat=0 do nlev=17,1,-1 ! from 8.16km until 0km p1 = pplay(i,nlev) if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then ! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 ! + ATB*epsilon50 + zeta50 ! Computation of the ATBperp of the phase discrimination line ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + & (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + & ATB(i,ncol,nlev)*epsilon50 + zeta50 !____________________________________________________________________________________________________ ! ! 4.2.a Ice: ATBperp above the phase discrimination line !____________________________________________________________________________________________________ ! ! ICE with temperature above 273,15°K = Liquid (false ice) if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds if(tmp(i,nlev).gt.273.15)then lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice ==> liq tmpl(i,ncol,nlev)=tmp(i,nlev) lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. cldlayphase(i,ncol,4,2) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,2) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,2) = 1. else ! low cloud cldlayphase(i,ncol,1,2) = 1. endif cldlayphase(i,ncol,4,5) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,5) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,5) = 1. else ! low cloud cldlayphase(i,ncol,1,5) = 1. endif else ! ICE with temperature below 273,15°K lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. tmpi(i,ncol,nlev)=tmp(i,nlev) cldlayphase(i,ncol,4,1) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,1) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,1) = 1. else ! low cloud cldlayphase(i,ncol,1,1) = 1. endif endif !____________________________________________________________________________________________________ ! ! 4.2.b Liquid: ATBperp below the phase discrimination line !____________________________________________________________________________________________________ ! else ! Liquid with temperature above 231,15°K if(tmp(i,nlev).gt.231.15)then lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. tmpl(i,ncol,nlev)=tmp(i,nlev) cldlayphase(i,ncol,4,2) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,2) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,2) = 1. else ! low cloud cldlayphase(i,ncol,1,2) = 1. endif else ! Liquid with temperature below 231,15°K = Ice (false liquid) tmpi(i,ncol,nlev)=tmp(i,nlev) lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liq ==> ice lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! false liq ==> ice cldlayphase(i,ncol,4,4) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,4) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,4) = 1. else ! low cloud cldlayphase(i,ncol,1,4) = 1. endif cldlayphase(i,ncol,4,1) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,1) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,1) = 1. else ! low cloud cldlayphase(i,ncol,1,1) = 1. endif endif endif ! end of discrimination condition toplvlsat=0 ! Find the level of the highest cloud with SR>30 if(x(i,ncol,nlev).gt.S_cld_att)then ! SR > 30. toplvlsat=nlev-1 goto 99 endif endif ! end of cloud condition enddo ! end of altitude loop 99 continue !____________________________________________________________________________________________________ ! ! Undefined phase: For a cloud located below another cloud with SR>30 ! see Cesana and Chepfer 2013 Sect.III.2 !____________________________________________________________________________________________________ ! if(toplvlsat.ne.0)then do nlev=toplvlsat,1,-1 p1 = pplay(i,nlev) if(cldy(i,ncol,nlev).eq.1.)then lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1. tmpu(i,ncol,nlev)=tmp(i,nlev) cldlayphase(i,ncol,4,3) = 1. ! tot cloud if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud cldlayphase(i,ncol,3,3) = 1. else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud cldlayphase(i,ncol,2,3) = 1. else ! low cloud cldlayphase(i,ncol,1,3) = 1. endif endif enddo endif toplvlsat=0 enddo enddo !____________________________________________________________________________________________________ ! ! Computation of final cloud phase diagnosis !____________________________________________________________________________________________________ ! ! Compute the Ice percentage in cloud = ice/(ice+liq) as a function ! of the occurrences lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2); WHERE (lidarcldphasetmp(:,:).gt. 0.) lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:) ELSEWHERE lidarcldphase(:,:,6) = undef ENDWHERE ! Compute Phase 3D Cloud Fraction WHERE ( nsub(:,:).gt.0.0 ) lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:) lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:) lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:) lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:) lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:) ELSEWHERE lidarcldphase(:,:,1) = undef lidarcldphase(:,:,2) = undef lidarcldphase(:,:,3) = undef lidarcldphase(:,:,4) = undef lidarcldphase(:,:,5) = undef ENDWHERE ! Compute Phase low mid high cloud fractions do iz = 1, Ncat do i=1,Nphase-3 do ic = 1, Ncolumns cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i) enddo enddo enddo do iz = 1, Ncat do i=4,5 do ic = 1, Ncolumns cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) enddo enddo enddo ! Compute the Ice percentage in cloud = ice/(ice+liq) cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2) WHERE (cldlayerphasetmp(:,:).gt. 0.) cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:) ELSEWHERE cldlayerphase(:,:,6) = undef ENDWHERE do i=1,Nphase-1 WHERE ( cldlayerphasesum(:,:).gt.0.0 ) cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:) ENDWHERE enddo do i=1,Npoints do iz=1,Ncat checkcldlayerphase=0. checkcldlayerphase2=0. if (cldlayerphasesum(i,iz).gt.0.0 )then do ic=1,Nphase-3 checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic) enddo checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz) endif enddo enddo do i=1,Nphase-1 WHERE ( nsublayer(:,:).eq.0.0 ) cldlayerphase(:,:,i) = undef ENDWHERE enddo ! Compute Phase 3D as a function of temperature do nlev=1,Nlevels do ncol=1,Ncolumns do i=1,Npoints do itemp=1,Ntemp if(tmpi(i,ncol,nlev).gt.0.)then if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1. endif elseif(tmpl(i,ncol,nlev).gt.0.)then if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1. endif elseif(tmpu(i,ncol,nlev).gt.0.)then if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1. endif endif enddo enddo enddo enddo ! Check temperature cloud fraction do i=1,Npoints do itemp=1,Ntemp checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4) if(checktemp.NE.lidarcldtemp(i,itemp,1))then print *, i,itemp print *, lidarcldtemp(i,itemp,1:4) endif enddo enddo ! Compute the Ice percentage in cloud = ice/(ice+liq) ! sumlidarcldtemp=sum(lidarcldtemp(:,:,2:3),3) sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3) WHERE(sumlidarcldtemp(:,:)>0.) lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:) ELSEWHERE lidarcldtemp(:,:,5)=undef ENDWHERE do i=1,4 WHERE(lidarcldtempind(:,:).gt.0.) lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:) ELSEWHERE lidarcldtemp(:,:,i) = undef ENDWHERE enddo RETURN END SUBROUTINE COSP_CLDFRAC ! --------------------------------------------------------------- ! BEGINNING OF OPAQ CHANGES ! #################################################################################### ! SUBROUTINE cosp_opaq ! Conventions: Ntype must be equal to 3 (opaque cloud, thin cloud, z_opaque) ! #################################################################################### SUBROUTINE COSP_OPAQ(Npoints,Ncolumns,Nlevels,Ntype,x,S_cld,undef,lidarcldtype, & cldtype,vgrid_z) IMPLICIT NONE ! Input arguments integer Npoints,Ncolumns,Nlevels,Ntype real x(Npoints,Ncolumns,Nlevels) real S_cld real undef real vgrid_z(Nlevels) ! Output : real lidarcldtype(Npoints,Nlevels,Ntype+1) ! 3D "lidar" OPAQ type + opacity fraction real cldtype(Npoints,Ntype) ! opaque and thin cloud covers, z_opaque ! Local variables integer ip, k, iz, ic, zopac real p1 real cldy(Npoints,Ncolumns,Nlevels) real cldyopaq(Npoints,Ncolumns,Nlevels) real srok(Npoints,Ncolumns,Nlevels) real srokopaq(Npoints,Ncolumns,Nlevels) real cldlay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover real nsublay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover real nsublayer(Npoints,Ntype) real nsub(Npoints,Nlevels) real nsubopaq(Npoints,Nlevels) real S_att_opaq real S_att ! #################################################################################### ! 1) Initialize ! #################################################################################### cldtype = 0.0 lidarcldtype = 0.0 nsub = 0.0 nsubopaq = 0.0 cldlay = 0.0 nsublay = 0.0 nsublayer = 0.0 S_att_opaq = 0.06 ! Fully Attenuated threshold, from Guzman et al. 2017, JGR-A S_att = 0.01 ! #################################################################################### ! 2) Cloud detection and Fully attenuated layer detection ! #################################################################################### do k=1,Nlevels ! Cloud detection at subgrid-scale: where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) ) cldy(:,:,k)=1.0 elsewhere cldy(:,:,k)=0.0 endwhere ! Fully attenuated layer detection at subgrid-scale: where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ne. undef) ) cldyopaq(:,:,k)=1.0 elsewhere cldyopaq(:,:,k)=0.0 endwhere ! Number of useful sub-column layers: where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) ) srok(:,:,k)=1.0 elsewhere srok(:,:,k)=0.0 endwhere ! Number of useful sub-columns layers for z_opaque 3D fraction: where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .ne. undef) ) srokopaq(:,:,k)=1.0 elsewhere srokopaq(:,:,k)=0.0 endwhere enddo ! #################################################################################### ! 3) Grid-box 3D OPAQ product fraction and cloud type cover (opaque/thin) + mean z_opaque ! #################################################################################### do k= Nlevels,1,-1 do ic = 1, Ncolumns do ip = 1, Npoints cldlay(ip,ic,1) = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k)) ! All clouds nsublay(ip,ic,1) = MAX(nsublay(ip,ic,1),srok(ip,ic,k)) nsublay(ip,ic,2) = MAX(nsublay(ip,ic,2),srok(ip,ic,k)) ! nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k)) nsub(ip,k) = nsub(ip,k) + srok(ip,ic,k) nsubopaq(ip,k) = nsubopaq(ip,k) + srokopaq(ip,ic,k) enddo enddo enddo ! OPAQ variables do ic = 1, Ncolumns do ip = 1, Npoints ! Declaring non-opaque cloudy profiles as thin cloud profiles if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then cldlay(ip,ic,2) = 1.0 endif ! Filling in 3D and 2D variables ! Opaque cloud profiles if ( cldlay(ip,ic,1) .eq. 1.0 ) then zopac = 0.0 do k=2,Nlevels ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0 cldlay(ip,ic,3) = vgrid_z(k-1) !z_opaque altitude nsublay(ip,ic,3) = 1.0 zopac = 1.0 endif if ( cldy(ip,ic,k) .eq. 1.0 ) then lidarcldtype(ip,k,1) = lidarcldtype(ip,k,1) + 1.0 endif enddo endif ! Thin cloud profiles if ( cldlay(ip,ic,2) .eq. 1.0 ) then do k=1,Nlevels ! Declaring thin cloud fraction for 3D variable if ( cldy(ip,ic,k) .eq. 1.0 ) then lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0 endif enddo endif enddo enddo ! 3D cloud types fraction (opaque=1 and thin=2) where ( nsub(:,:) .gt. 0.0 ) lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:) lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:) elsewhere lidarcldtype(:,:,1) = undef lidarcldtype(:,:,2) = undef endwhere ! 3D z_opaque fraction (=3) where ( nsubopaq(:,:) .gt. 0.0 ) lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:) elsewhere lidarcldtype(:,:,3) = undef endwhere ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1) lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3) do ip = 1, Npoints do k = Nlevels-1, 1, -1 if ( lidarcldtype(ip,k,3) .ne. undef ) then lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3) endif enddo enddo where ( nsubopaq(:,:) .eq. 0.0 ) lidarcldtype(:,:,4) = undef endwhere ! Layered cloud types (opaque, thin and z_opaque 2D variables) do iz = 1, Ntype do ic = 1, Ncolumns cldtype(:,iz) = cldtype(:,iz) + cldlay(:,ic,iz) nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz) enddo enddo where (nsublayer(:,:) .gt. 0.0) cldtype(:,:) = cldtype(:,:)/nsublayer(:,:) elsewhere cldtype(:,:) = undef endwhere END SUBROUTINE COSP_OPAQ ! END OF OPAQ CHANGES END MODULE MOD_LMD_IPSL_STATS