Changeset 2435 for LMDZ5/branches/testing/libf/phylmd/cosp/lidar_simulator.F90

Timestamp:

Jan 28, 2016, 5:02:13 PM (9 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2396:2434 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/cosp/lidar_simulator.F90 (modified) (18 diffs)

LMDZ5/branches/testing/libf/phylmd/cosp/lidar_simulator.F90

@@ -1 +1 @@
 ! 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/lidar_simulator.F90 $
 ! 
 ! Redistribution and use in source and binary forms, with or without modification, are permitted 
@@ -28 +30 @@
                 , q_lsliq, q_lsice, q_cvliq, q_cvice &
                 , ls_radliq, ls_radice, cv_radliq, cv_radice &
-                , frac_out, ice_type &
-                , pmol, pnorm, tautot, refl )
+                , ice_type, pmol, pnorm, pnorm_perp_tot,tautot, refl )
 !
 !---------------------------------------------------------------------------------
@@ -75 +76 @@
 ! - Bug fix in computation of pmol and pnorm, thanks to Masaki Satoh: a factor 2 
 ! was missing. This affects the ATB values but not the cloud fraction. 
+!
+! January 2013, G. Cesana and H. Chepfer:
+! - Add the perpendicular component of the backscattered signal (pnorm_perp_tot) in the arguments
+! - Add the temperature for each levels (temp) in the arguments
+! - Add the computation of the perpendicular component of the backscattered lidar signal 
+! 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
 !
 !---------------------------------------------------------------------------------
@@ -96 +104 @@
 !  cv_radliq: effective radius of CONV liquid particles (meters)
 !  cv_radice: effective radius of CONV ice particles (meters)
-!  frac_out : cloud cover in each sub-column of the gridbox (output from scops)
 !  ice_type : ice particle shape hypothesis (ice_type=0 for spheres, ice_type=1 
 !             for non spherical particles)
@@ -103 +110 @@
 !  pmol : molecular attenuated backscatter lidar signal power (m^-1.sr^-1)
 !  pnorm: total attenuated backscatter lidar signal power (m^-1.sr^-1)
+!  pnorm_perp_tot: perpendicular attenuated backscatter lidar signal power (m^-1.sr^-1)
 !  tautot: optical thickess integrated from top to level z
 !  refl : parasol(polder) reflectance
@@ -130 +138 @@
       REAL pres(npoints,nlev)    ! pressure full levels
       REAL presf(npoints,nlev+1) ! pressure half levels
-      REAL temp(npoints,nlev)
       REAL q_lsliq(npoints,nlev), q_lsice(npoints,nlev)
       REAL q_cvliq(npoints,nlev), q_cvice(npoints,nlev)
       REAL ls_radliq(npoints,nlev), ls_radice(npoints,nlev)
       REAL cv_radliq(npoints,nlev), cv_radice(npoints,nlev)
-      REAL frac_out(npoints,nlev)
 
 ! outputs (for each subcolumn):
@@ -168 +174 @@
 
 !   sub-column variables:
-      REAL frac_sub(npoints,nlev)
       REAL qpart(npoints,nlev,npart) ! mixing ratio particles in each subcolumn
       REAL alpha_part(npoints,nlev,npart)
@@ -177 +182 @@
       REAL tautot_lay(npoints)   ! temporary variable, total opt. thickness of layer k
 !     Optical thickness from TOA to surface for Parasol
-      REAL tautot_S_liq(npoints),tautot_S_ice(npoints)     ! for liq and ice clouds
-
-! Abderrahmane 8-2-2011
-      Logical iflag_testlidar
-      PARAMETER (iflag_testlidar=.false.)
+     REAL tautot_S_liq(npoints),tautot_S_ice(npoints)     ! for liq and ice clouds
+
+
+! Local variables
+      REAL Alpha, Beta, Gamma  ! Polynomial coefficient for ATBperp computation
+      REAL temp(npoints,nlev)                   ! temperature of layer k
+      REAL betatot_ice(npoints,nlev)    ! backscatter coefficient for ice particles
+      REAL beta_perp_ice(npoints,nlev)  ! perpendicular backscatter coefficient for ice
+      REAL betatot_liq(npoints,nlev)    ! backscatter coefficient for liquid particles
+      REAL beta_perp_liq(npoints,nlev)  ! perpendicular backscatter coefficient for liq
+      REAL tautot_ice(npoints,nlev)     ! total optical thickness of ice
+      REAL tautot_liq(npoints,nlev)     ! total optical thickness of liq
+      REAL tautot_lay_ice(npoints)    ! total optical thickness of ice in the layer k
+      REAL tautot_lay_liq(npoints)    ! total optical thickness of liq in the layer k
+      REAL pnorm_liq(npoints,nlev)    ! lidar backscattered signal power for liquid
+      REAL pnorm_ice(npoints,nlev)    ! lidar backscattered signal power for ice
+      REAL pnorm_perp_ice(npoints,nlev) ! perpendicular lidar backscattered signal power for ice
+      REAL pnorm_perp_liq(npoints,nlev) ! perpendicular lidar backscattered signal power for liq
+
+! Output variable
+      REAL pnorm_perp_tot (npoints,nlev) ! perpendicular lidar backscattered signal power
+
+!------------------------------------------------------------
+!---- 0. Initialisation :
+!------------------------------------------------------------
+betatot_ice(:,:)=0
+betatot_liq(:,:)=0
+beta_perp_ice(:,:)=0
+beta_perp_liq(:,:)=0
+tautot_ice(:,:)=0
+tautot_liq(:,:)=0
+tautot_lay_ice(:)=0;
+tautot_lay_liq(:)=0;
+pnorm_liq(:,:)=0
+pnorm_ice(:,:)=0
+pnorm_perp_ice(:,:)=0
+pnorm_perp_liq(:,:)=0
+pnorm_perp_tot(:,:)=0
+
+
+! Polynomial coefficients (Alpha, Beta, Gamma) which allow to compute the ATBperpendicular
+! as a function of the ATB for ice or liquid cloud particles derived from CALIPSO-GOCCP
+! observations at 120m vertical grid (Cesana and Chepfer, JGR, 2013).
+!
+! Relationship between ATBice and ATBperp,ice for ice particles
+!  ATBperp,ice = Alpha*ATBice 
+         Alpha = 0.2904
+
+! Relationship between ATBice and ATBperp,ice for liquid particles
+!  ATBperp,ice = Beta*ATBice^2 + Gamma*ATBice
+         Beta = 0.4099
+         Gamma = 0.009
 
 !------------------------------------------------------------
@@ -201 +253 @@
 ! We repeat the same coefficients for LS and CONV cloud to make code more readable
 !*     LS Liquid water coefficients:
-         polpart(INDX_LSLIQ,1) =  2.6980e-8     
+         polpart(INDX_LSLIQ,1) =  2.6980e-8
          polpart(INDX_LSLIQ,2) = -3.7701e-6
          polpart(INDX_LSLIQ,3) =  1.6594e-4
@@ -208 +260 @@
 !*     LS Ice coefficients: 
       if (ice_type.eq.0) then     
-         polpart(INDX_LSICE,1) = -1.0176e-8   
+         polpart(INDX_LSICE,1) = -1.0176e-8
          polpart(INDX_LSICE,2) =  1.7615e-6
          polpart(INDX_LSICE,3) = -1.0480e-4
@@ -216 +268 @@
 !*     LS Ice NS coefficients: 
       if (ice_type.eq.1) then 
-         polpart(INDX_LSICE,1) = 1.3615e-8  
-         polpart(INDX_LSICE,2) = -2.04206e-6 
+         polpart(INDX_LSICE,1) = 1.3615e-8
+         polpart(INDX_LSICE,2) = -2.04206e-6
          polpart(INDX_LSICE,3) = 7.51799e-5
          polpart(INDX_LSICE,4) = 0.00078213
@@ -223 +275 @@
       endif
 !*     CONV Liquid water coefficients:
-         polpart(INDX_CVLIQ,1) =  2.6980e-8     
+         polpart(INDX_CVLIQ,1) =  2.6980e-8
          polpart(INDX_CVLIQ,2) = -3.7701e-6
          polpart(INDX_CVLIQ,3) =  1.6594e-4
@@ -230 +282 @@
 !*     CONV Ice coefficients: 
       if (ice_type.eq.0) then 
-         polpart(INDX_CVICE,1) = -1.0176e-8   
+         polpart(INDX_CVICE,1) = -1.0176e-8
          polpart(INDX_CVICE,2) =  1.7615e-6
          polpart(INDX_CVICE,3) = -1.0480e-4
@@ -268 +320 @@
                   -(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
       enddo
-
-! cloud fraction (0 or 1) in each sub-column:
-! (if frac_out=1or2 -> frac_sub=1; if frac_out=0 -> frac_sub=0)
-      frac_sub = MIN( frac_out, 1.0 )
 
 !------------------------------------------------------------
@@ -316 +364 @@
 
 !------------------------------------------------------------
-!---- 4. Backscatter signal:
+!---- 4.1 Total Backscatter signal:
 !------------------------------------------------------------
 
@@ -356 +404 @@
         END WHERE
       END DO
-!
+
 ! Total signal (molecular + particules):
+!
 !
 ! For performance reason on vector computers, the 2 following lines should not be used
@@ -373 +422 @@
       pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
             & * (1.-exp(-2.0*tautot(:,nlev)))
+
 !     Other layers
       DO k= nlev-1, 1, -1
-        tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
+          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
         WHERE (tautot_lay(:).GT.0.)
-       pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
-!correc          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) & ! correc Satoh
-!correc               &               / (2.0*tautot_lay(:)) &          ! correc Satoh
+          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
                & * (1.-EXP(-2.0*tautot_lay(:)))
         ELSEWHERE
@@ -387 +435 @@
       END DO
 
-     if (iflag_testlidar) then
-!+JLD test 
-!     do k=1,nlev
-!      print*,'Min val de frac_out=',k,minval(frac_out(:,k))
-!      print*,'Max val de frac_out=',k,maxval(frac_out(:,k))
-!     enddo
-       where ( frac_out(:,:).ge.0.5)
-! Correction AI 9 5 11          pnorm(:,:) = pmol(:,:)*10.
-       pnorm(:,:) = pmol(:,:)*50.
-        elsewhere
-          pnorm(:,:) = pmol(:,:)
-        endwhere
-!-JLD test 
-     endif
+!------------------------------------------------------------
+!---- 4.2 Ice/Liq Backscatter signal:
+!------------------------------------------------------------
+
+! Contribution of the molecular to beta
+      betatot_ice(:,:) = beta_mol(:,:)
+      betatot_liq(:,:) = beta_mol(:,:)
+
+      tautot_ice(:,:) = tau_mol(:,:)
+      tautot_liq(:,:) = tau_mol(:,:)
+
+      DO i = 2, npart,2
+           betatot_ice(:,:) = betatot_ice(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
+           tautot_ice(:,:) = tautot_ice(:,:)  + tau_part(:,:,i)
+      ENDDO ! i
+      DO i = 1, npart,2
+           betatot_liq(:,:) = betatot_liq(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
+           tautot_liq(:,:) = tautot_liq(:,:)  + tau_part(:,:,i)
+      ENDDO ! i
+
+
+! Computation of the ice and liquid lidar backscattered signal (ATBice and ATBliq)
+!     Ice only
+!     Upper layer
+      pnorm_ice(:,nlev) = betatot_ice(:,nlev) / (2.*tautot_ice(:,nlev)) &
+            & * (1.-exp(-2.0*tautot_ice(:,nlev)))
+
+      DO k= nlev-1, 1, -1
+          tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1) 
+        WHERE (tautot_lay_ice(:).GT.0.)
+         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))/(2.*tautot_lay_ice(:)) &
+               & * (1.-EXP(-2.0*tautot_lay_ice(:)))
+        ELSEWHERE
+         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))
+        END WHERE
+      ENDDO
+
+!     Liquid only
+!     Upper layer
+      pnorm_liq(:,nlev) = betatot_liq(:,nlev) / (2.*tautot_liq(:,nlev)) &
+            & * (1.-exp(-2.0*tautot_liq(:,nlev)))
+
+      DO k= nlev-1, 1, -1
+          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1) 
+        WHERE (tautot_lay_liq(:).GT.0.)
+          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))/(2.*tautot_lay_liq(:)) &
+               & * (1.-EXP(-2.0*tautot_lay_liq(:)))
+        ELSEWHERE
+          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))
+        END WHERE
+      ENDDO
+
+
+! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering 
+! contribution (Cesana and Chepfer 2013, JGR)
+!  ATBperp,ice = Alpha*ATBice 
+!  ATBperp,liq = Beta*ATBliq^2 + Gamma*ATBliq
+
+      DO k= nlev, 1, -1
+              pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
+              pnorm_perp_liq(:,k) = 1000*Beta * pnorm_liq(:,k)**2 + Gamma * pnorm_liq(:,k) ! Liquid particles
+      ENDDO
+
+! Computation of beta_perp_ice/liq using the lidar equation
+!     Ice only
+!     Upper layer 
+      beta_perp_ice(:,nlev) = pnorm_perp_ice(:,nlev) * (2.*tautot_ice(:,nlev)) &
+            & / (1.-exp(-2.0*tautot_ice(:,nlev)))
+
+      DO k= nlev-1, 1, -1
+        tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
+        WHERE (tautot_lay_ice(:).GT.0.)
+         beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
+            & / (1.-exp(-2.0*tautot_lay_ice(:)))
+
+        ELSEWHERE
+         beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/EXP(-2.0*tautot_ice(:,k+1))
+        END WHERE
+      ENDDO
+
+!     Liquid only
+!     Upper layer 
+      beta_perp_liq(:,nlev) = pnorm_perp_liq(:,nlev) * (2.*tautot_liq(:,nlev)) &
+            & / (1.-exp(-2.0*tautot_liq(:,nlev)))
+
+      DO k= nlev-1, 1, -1
+          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1) 
+        WHERE (tautot_lay_liq(:).GT.0.)
+         beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ EXP(-2.0*tautot_liq(:,k+1)) * (2.*tautot_lay_liq(:)) &
+            & / (1.-exp(-2.0*tautot_lay_liq(:)))
+
+        ELSEWHERE
+         beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/EXP(-2.0*tautot_liq(:,k+1))
+        END WHERE
+      ENDDO
+
+
+
+!------------------------------------------------------------
+!---- 4.3 Perpendicular Backscatter signal:
+!------------------------------------------------------------
+
+! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
+!     Upper layer 
+    WHERE(tautot(:,nlev).GT.0)
+          pnorm_perp_tot(:,nlev) = &
+              (beta_perp_ice(:,nlev)+beta_perp_liq(:,nlev)-(beta_mol(:,nlev)/(1+1/0.0284))) / (2.*tautot(:,nlev)) &
+              & * (1.-exp(-2.0*tautot(:,nlev)))
+    ELSEWHERE
+    pnorm_perp_tot(:,nlev) = 0.
+    ENDWHERE
+
+!     Other layers
+      DO k= nlev-1, 1, -1
+          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
+
+          ! The perpendicular component of the molecular backscattered signal (Betaperp) has been 
+          ! taken into account two times (once for liquid and once for ice). 
+          ! We remove one contribution using 
+          ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following equations:
+            WHERE (pnorm(:,k).eq.0)
+                  pnorm_perp_tot(:,k)=0.
+                  ELSEWHERE
+                    WHERE (tautot_lay(:).GT.0.)
+                      pnorm_perp_tot(:,k) = &
+                          (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
+                          EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
+                          & * (1.-EXP(-2.0*tautot_lay(:)))
+                    ELSEWHERE
+          !         This must never happen, but just in case, to avoid div. by 0
+                      pnorm_perp_tot(:,k) = &
+                           (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
+                          EXP(-2.0*tautot(:,k+1))
+                    END WHERE
+            ENDWHERE
+
+      END DO
 
 !-------- End computation Lidar --------------------------