Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90	(revision 5158)
@@ -78,9 +78,9 @@
     box_MISR_ztop(1:npoints,1:ncol) = 0._wp  
 
-    do j=1,npoints
+    DO j=1,npoints
 
        ! Estimate distribution of Model layer tops
        dist_model_layertops(j,:)=0
-       do ilev=1,nlev
+       DO ilev=1,nlev
           ! Define location of "layer top"
           if(ilev.eq.1 .or. ilev.eq.nlev) then
@@ -93,5 +93,5 @@
           ! *NOTE* the first MISR level is "no height" level
           iMISR_ztop=2
-          do loop=2,numMISRHgtBins
+          DO loop=2,numMISRHgtBins
              if ( ztest .gt. 1000*misr_histHgt(loop+1) ) then
                 iMISR_ztop=loop+1
@@ -103,10 +103,10 @@
 
        ! For each GCM cell or horizontal model grid point   
-       do ibox=1,ncol
+       DO ibox=1,ncol
           ! Compute optical depth as a cummulative distribution in the vertical (nlev).
           tauOUT(j,ibox)=sum(dtau(j,ibox,1:nlev))
 
           thres_crossed_MISR=0
-          do ilev=1,nlev
+          DO ilev=1,nlev
              ! If there a cloud, start the counter and store this height
              if(thres_crossed_MISR .eq. 0 .and. dtau(j,ibox,ilev) .gt. 0.) then
@@ -187,5 +187,5 @@
     ! This setup assumes the columns represent a about a 1 to 4 km scale
     ! it will need to be modified significantly, otherwise
-!	! DS2015: Add loop over gridpoints and index accordingly.
+!    ! DS2015: Add loop over gridpoints and index accordingly.
 !    if(ncol.eq.1) then
 !       ! Adjust based on neightboring points.
@@ -214,5 +214,5 @@
      
     ! Fill dark scenes 
-    do j=1,numMISRHgtBins
+    DO j=1,numMISRHgtBins
        where(sunlit .ne. 1) dist_model_layertops(1:npoints,j) = R_UNDEF
     enddo
@@ -254,5 +254,5 @@
     tauWRK(1:npoints,1:ncol)                      = tau(1:npoints,1:ncol)
     box_MISR_ztopWRK(1:npoints,1:ncol)            = box_MISR_ztop(1:npoints,1:ncol)
-    do j=1,npoints
+    DO j=1,npoints
 
        ! Subcolumns that are cloudy(true) and not(false)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90	(revision 5158)
@@ -239,5 +239,5 @@
      real(wp),dimension(:),pointer :: &
           isccp_totalcldarea => null(), & ! The fraction of model grid box columns with cloud
-           		                  ! somewhere in them. (%)
+                                     ! somewhere in them. (%)
           isccp_meantb => null(),       & ! Mean all-sky 10.5 micron brightness temperature. (K)
           isccp_meantbclr => null(),    & ! Mean clear-sky 10.5 micron brightness temperature. (K)
@@ -837,5 +837,5 @@
        allocate(parasolPix_refl(parasolIN%Npoints,parasolIN%Ncolumns,PARASOL_NREFL))
        ! Call simulator
-       do icol=1,parasolIN%Ncolumns
+       DO icol=1,parasolIN%Ncolumns
           call parasol_subcolumn(parasolIN%npoints, PARASOL_NREFL,                       &
                                  parasolIN%tautot_S_liq(1:parasolIN%Npoints,icol),       &
@@ -855,5 +855,5 @@
        allocate(cloudsatDBZe(cloudsatIN%Npoints,cloudsatIN%Ncolumns,cloudsatIN%Nlevels), &
                 cloudsatZe_non(cloudsatIN%Npoints,cloudsatIN%Ncolumns,cloudsatIN%Nlevels))
-       do icol=1,cloudsatIN%ncolumns
+       DO icol=1,cloudsatIN%ncolumns
           call quickbeam_subcolumn(cloudsatIN%rcfg,cloudsatIN%Npoints,cloudsatIN%Nlevels,&
                                    cloudsatIN%hgt_matrix/1000._wp,                       &
@@ -876,5 +876,5 @@
                    modisRetrievedCloudTopPressure(modisIN%nSunlit,modisIN%nColumns))
           ! Call simulator
-          do i = 1, modisIN%nSunlit
+          DO i = 1, modisIN%nSunlit
              call modis_subcolumn(modisIN%Ncolumns,modisIN%Nlevels,modisIN%pres(i,:),    &
                                   modisIN%tau(int(modisIN%sunlit(i)),:,:),               &
@@ -1246,5 +1246,5 @@
 !            cospgridIN%hgt_matrix, cospgridIN%hgt_matrix_half,                           & !PREC_BUG
             cospgridIN%land(:), cospgridIN%surfelev(:), cospgridIN%at(:,cospIN%Nlevels), & !PREC_BUG
-	    cospIN%fracPrecipIce, cospgridIN%hgt_matrix, cospgridIN%hgt_matrix_half,     & !PREC_BUG
+        cospIN%fracPrecipIce, cospgridIN%hgt_matrix, cospgridIN%hgt_matrix_half,     & !PREC_BUG
             cospOUT%cloudsat_cfad_ze(ij:ik,:,:), cospOUT%cloudsat_precip_cover,          &
             cospOUT%cloudsat_pia)
@@ -1691,5 +1691,5 @@
        ! Other grid requested. Constant vertical spacing with top at 20 km
        if (.not. luseCSATvgrid) zstep = 20000._wp/Nvgrid
-       do i=1,Nvgrid
+       DO i=1,Nvgrid
           vgrid_zl(Nlvgrid-i+1) = (i-1)*zstep
           vgrid_zu(Nlvgrid-i+1) = i*zstep
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_atlid_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_atlid_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_atlid_interface.F90	(revision 5158)
@@ -62,5 +62,5 @@
 
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !	END MODULE
+  !    END MODULE
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 END MODULE MOD_COSP_ATLID_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_calipso_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_calipso_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_calipso_interface.F90	(revision 5158)
@@ -81,5 +81,5 @@
 
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !	END MODULE
+  !    END MODULE
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 END MODULE MOD_COSP_CALIPSO_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_cloudsat_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_cloudsat_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_cloudsat_interface.F90	(revision 5158)
@@ -118,5 +118,5 @@
     ! Set up Re bin "structure" for z_scaling
     rcfg%base_list(1)=0
-    do j=1,Re_MAX_BIN
+    DO j=1,Re_MAX_BIN
        rcfg%step_list(j)=0.1_wp+0.1_wp*((j-1)**1.5)
        if(rcfg%step_list(j)>Re_BIN_LENGTH) then
@@ -138,5 +138,5 @@
 
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !                              	  END MODULE
+  !                                    END MODULE
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 END MODULE MOD_COSP_CLOUDSAT_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_config.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_config.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_config.F90	(revision 5158)
@@ -180,5 +180,5 @@
                                    1.0,      0.909013, 0.709554, 0.430405, 0.121567/),   &
                                    shape=(/PARASOL_NREFL,PARASOL_NTAU/)),                & 
-         ! LUT for ice particles         			     
+         ! LUT for ice particles
          rlumB = reshape(source=(/ 0.03,     0.03,     0.03,     0.03,     0.03,         &
                                    0.092170, 0.087082, 0.083325, 0.084935, 0.054157,     &
@@ -422,8 +422,8 @@
          vgrid_z_in = (/240.0, 720.0, 1200.0, 1680.0, 2160.0, 2640.0, 3120.0, 3600.0, &
                         4080.0, 4560.0, 5040.0, 5520.0, 6000.0, 6480.0, 6960.0, 7440.0, &
-			7920.0, 8400.0, 8880.0, 9360.0, 9840.0, 10320.0, 10800.0, &
-			11280.0, 11760.0, 12240.0, 12720.0, 13200.0, 13680.0, 14160.0, &
-			14640.0, 15120.0, 15600.0, 16080.0, 16560.0, 17040.0, 17520.0, &
-			18000.0, 18480.0, 18960.0/)
+            7920.0, 8400.0, 8880.0, 9360.0, 9840.0, 10320.0, 10800.0, &
+            11280.0, 11760.0, 12240.0, 12720.0, 13200.0, 13680.0, 14160.0, &
+            14640.0, 15120.0, 15600.0, 16080.0, 16560.0, 17040.0, 17520.0, &
+            18000.0, 18480.0, 18960.0/)
 
 END MODULE MOD_COSP_CONFIG
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_grLidar532_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_grLidar532_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_grLidar532_interface.F90	(revision 5158)
@@ -58,5 +58,5 @@
   
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !	END MODULE
+  !    END MODULE
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 END MODULE MOD_COSP_GRLIDAR532_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_isccp_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_isccp_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_isccp_interface.F90	(revision 5158)
@@ -67,5 +67,5 @@
 
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !  							SUBROUTINE cosp_isccp_init
+  !                              SUBROUTINE cosp_isccp_init
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   SUBROUTINE COSP_ISCCP_INIT(top_height,top_height_direction)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_misr_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_misr_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_misr_interface.F90	(revision 5158)
@@ -36,5 +36,5 @@
 
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  ! 								TYPE misr_in
+  !                                 TYPE misr_in
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   type misr_IN
@@ -56,5 +56,5 @@
 
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !  							SUBROUTINE cosp_misr_init
+  !                              SUBROUTINE cosp_misr_init
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   SUBROUTINE COSP_MISR_INIT()
@@ -63,5 +63,5 @@
 
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  ! 							    	END MODULE
+  !                                     END MODULE
   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 END MODULE MOD_COSP_MISR_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90	(revision 5158)
@@ -71,5 +71,5 @@
     
     varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
-    do j=1,dim2
+    DO j=1,dim2
        where(flag(:,j,:) .eq. 1)
           varOUT(:,j,:) = varIN2
@@ -135,5 +135,5 @@
     
     
-    do i=1,npoints
+    DO i=1,npoints
        where(cloudIce(i,:, :) <= 0.) 
           fracL(:, :) = 1._wp
@@ -170,6 +170,6 @@
     w0(1:nPoints,1:nSubCols,1:nLevels) = 0._wp
     
-    do j =1,nPoints
-       do i=1,nSubCols
+    DO j =1,nPoints
+       DO i=1,nSubCols
           water_g(1:nLevels)  = get_g_nir(  phaseIsLiquid, sizeLIQ(j,i,1:nLevels)) 
           water_w0(1:nLevels) = get_ssa_nir(phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
@@ -189,5 +189,5 @@
     
     ! Compute the total optical thickness and the proportion due to liquid in each cell
-    do i=1,npoints
+    DO i=1,npoints
        where(tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels) > 0.) 
           fracLIQ(i,1:nSubCols,1:nLevels) = tauLIQ(i,1:nSubCols,1:nLevels)/ &
@@ -358,5 +358,5 @@
     ! Altitude at half pressure levels:
     zheight(1:npoints,nlev+1) = 0._wp
-    do k=nlev,1,-1
+    DO k=nlev,1,-1
        zheight(1:npoints,k) = zheight(1:npoints,k+1) &
             -(presf(1:npoints,k)-presf(1:npoints,k+1))/(rhoair(1:npoints,k)*grav)
@@ -392,5 +392,5 @@
     ! ##############################################################################
     ! Polynomials kp_lidar derived from Mie theory
-    do i = 1, npart
+    DO i = 1, npart
        where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
           kp_part(1:npoints,1:nlev,i) = &
@@ -412,5 +412,5 @@
 
     ! Loop over all subcolumns
-    do icol=1,ncolumns
+    DO icol=1,ncolumns
        ! ##############################################################################
        ! Mixing ratio particles in each subcolum
@@ -425,5 +425,5 @@
        ! ##############################################################################
        ! Alpha of particles in each subcolumn:
-       do i = 1, npart
+       DO i = 1, npart
           where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
              alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat &
@@ -437,10 +437,10 @@
        ! Optical thicknes
        tau_part(1:npoints,1:nlev,1:npart) = rdiffm * alpha_part(1:npoints,1:nlev,1:npart)
-       do i = 1, npart
+       DO i = 1, npart
           ! Optical thickness of each layer (particles)
           tau_part(1:npoints,1:nlev,i) = tau_part(1:npoints,1:nlev,i) &
                & * (zheight(1:npoints,1:nlev)-zheight(1:npoints,2:nlev+1) )
           ! Optical thickness from TOA to layer k (particles)
-          do k=zi,zf,zinc
+          DO k=zi,zf,zinc
              tau_part(1:npoints,k,i) = tau_part(1:npoints,k,i) + tau_part(1:npoints,k+zoffset,i)
           enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_parasol_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_parasol_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_parasol_interface.F90	(revision 5158)
@@ -35,5 +35,5 @@
 
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !									TYPE cosp_parasol
+  !                                    TYPE cosp_parasol
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   TYPE PARASOL_SGX
@@ -63,5 +63,5 @@
   END TYPE COSP_PARASOL
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  !										TYPE parasol_in
+  !                                        TYPE parasol_in
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   TYPE parasol_IN
@@ -86,5 +86,5 @@
 
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  ! 								    END MODULE
+  !                                     END MODULE
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 end module MOD_COSP_PARASOL_INTERFACE
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90	(revision 5158)
@@ -81,5 +81,5 @@
    r = 0._wp
 
-   do i=1,Npoints
+   DO i=1,Npoints
      ! Calculate tops and bottoms of new and old grids
      oldgrid_bot = zhalf(i,:)
@@ -88,9 +88,9 @@
      l = 0 ! Index of level in the old grid
      ! Loop over levels in the new grid
-     do k = 1,Nglevels
+     DO k = 1,Nglevels
        Nw = 0 ! Number of weigths
        wt = 0._wp ! Sum of weights
        ! Loop over levels in the old grid and accumulate total for weighted average
-       do
+       DO
          l = l + 1
          w = 0.0 ! Initialise weight to 0
@@ -119,5 +119,5 @@
              Nw = Nw + 1
              wt = wt + w
-             do j=1,Ncolumns
+             DO j=1,Ncolumns
                if (lunits) then
                  if (y(i,j,l) /= R_UNDEF) then
@@ -138,5 +138,5 @@
        ! Calculate average in new grid
        if (Nw > 0) then
-         do j=1,Ncolumns
+         DO j=1,Ncolumns
            r(i,j,k) = r(i,j,k)/wt
          enddo
@@ -146,7 +146,7 @@
 
    ! Set points under surface to R_UNDEF, and change to dBZ if necessary
-   do k=1,Nglevels
-     do j=1,Ncolumns
-       do i=1,Npoints
+   DO k=1,Nglevels
+     DO j=1,Ncolumns
+       DO i=1,Npoints
          if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
            if (lunits) then
@@ -202,6 +202,6 @@
     radar_tcc = 0._wp
     radar_tcc2 = 0._wp
-    do pr=1,Npoints
-       do i=1,Ncolumns
+    DO pr=1,Npoints
+       DO i=1,Ncolumns
           flag_sat = 0
           flag_cld = 0
@@ -210,10 +210,10 @@
           ! look for j_1km from bottom to top
           j = 1
-          do while (Ze_tot(pr,i,j) .eq. R_GROUND)
+          DO while (Ze_tot(pr,i,j) .eq. R_GROUND)
              j = j+1
           enddo
           j_1km = j+1  !this is the vertical index of 1km above surface  
           
-          do j=1,Nlevels
+          DO j=1,Nlevels
              sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
              if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
@@ -266,5 +266,5 @@
     integer :: ij
     
-    do ij=2,Nbins+1  
+    DO ij=2,Nbins+1
        hist1D(ij-1) = count(var .ge. bins(ij-1) .and. var .lt. bins(ij))
        if (count(var .eq. R_GROUND) .ge. 1) hist1D(ij-1)=R_UNDEF
@@ -298,6 +298,6 @@
     integer :: ij,ik
     
-    do ij=2,nbin1+1
-       do ik=2,nbin2+1
+    DO ij=2,nbin1+1
+       DO ik=2,nbin2+1
           jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
                var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))        
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_utils.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_utils.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_utils.F90	(revision 5158)
@@ -77,7 +77,7 @@
         delta = (alpha_x + b_x + d_x - n_bx + 1._wp)
         
-        do k=1,Nlevels
-            do j=1,Ncolumns
-                do i=1,Npoints
+        DO k=1,Nlevels
+            DO j=1,Ncolumns
+                DO i=1,Npoints
                     if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then
                         rho = p(i,k)/(287.05_wp*T(i,k))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90	(revision 5158)
@@ -135,8 +135,8 @@
     
     if (debugcol.ne.0) then
-       do j=1,npoints,debugcol
+       DO j=1,npoints,debugcol
           
           ! Produce character output
-          do ilev=1,nlev
+          DO ilev=1,nlev
              acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2
              where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
@@ -151,5 +151,5 @@
           write(9,'(a1)') ' '
           
-          do ibox=1,ncol
+          DO ibox=1,ncol
              write(9,'(40(a1),1x,40(a1))') &
                   (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),&
@@ -231,5 +231,5 @@
        itrop(1:npoints)     = 1
 
-       do ilev=1,nlev
+       DO ilev=1,nlev
           where(pfull(1:npoints,ilev) .lt. 40000. .and. &
                 pfull(1:npoints,ilev) .gt.  5000. .and. &
@@ -242,5 +242,5 @@
        enddo
 
-       do ilev=1,nlev
+       DO ilev=1,nlev
           atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),&
                at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev  .ge. itrop(1:npoints))
@@ -256,5 +256,5 @@
        ! at a wavenumber of 955 cm-1, or 10.47 microns 
        ! ############################################################################
-       do ilev=1,nlev
+       DO ilev=1,nlev
           press(1:npoints)  = pfull(1:npoints,ilev)*10._wp
           dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10
@@ -273,5 +273,5 @@
        fluxtop_clrsky(1:npoints)            = 0._wp
        trans_layers_above_clrsky(1:npoints) = 1._wp
-       do ilev=1,nlev
+       DO ilev=1,nlev
           ! Black body emission at temperature of the layer
           bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp )
@@ -300,9 +300,9 @@
        fluxtop(1:npoints,1:ncol)            = 0._wp
        trans_layers_above(1:npoints,1:ncol) = 1._wp
-       do ilev=1,nlev
+       DO ilev=1,nlev
           ! Black body emission at temperature of the layer
           bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp)
           
-          do ibox=1,ncol
+          DO ibox=1,ncol
              ! Emissivity
              dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
@@ -320,5 +320,5 @@
        ! Add in surface emission
        bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
-       do ibox=1,ncol
+       DO ibox=1,ncol
           fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lw*bb(1:npoints)*trans_layers_above(1:npoints,ibox) 
        end do
@@ -344,10 +344,10 @@
        btcmin(1:npoints) = 1._wp /  ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp ) 
 
-       do ibox=1,ncol
+       DO ibox=1,ncol
           transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints))
           tauir(1:npoints)    = tau(1:npoints,ibox)/2.13_wp
           taumin(1:npoints)   = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
           if (isccp_top_height .eq. 1) then
-             do j=1,npoints  
+             DO j=1,npoints
                 if (transmax(j) .gt. 0.001 .and.  transmax(j) .le. 0.9999999) then
                    fluxtopinit(j) = fluxtop(j,ibox)
@@ -355,6 +355,6 @@
                 endif
              enddo
-             do icycle=1,2
-                do j=1,npoints  
+             DO icycle=1,2
+                DO j=1,npoints
                    if (tau(j,ibox) .gt. (tauchk)) then 
                       if (transmax(j) .gt. 0.001 .and.  transmax(j) .le. 0.9999999) then
@@ -397,5 +397,5 @@
     ! pressure (isccp_top_height = 1 or 3)
     ! ####################################################################################
-    do ibox=1,ncol
+    DO ibox=1,ncol
        !segregate according to optical thickness
        if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then  
@@ -403,7 +403,7 @@
           ! Find level whose temperature most closely matches brightness temperature
           nmatch(1:npoints)=0
-          do k1=1,nlev-1
+          DO k1=1,nlev-1
              ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)        
-             do j=1,npoints 
+             DO j=1,npoints
                 if (ilev           .ge. itrop(j)     .and. &
                      ((at(j,ilev)  .ge. tb(j,ibox)   .and. &  
@@ -417,5 +417,5 @@
           enddo
 
-          do j=1,npoints 
+          DO j=1,npoints
              if (nmatch(j) .ge. 1) then
                 k1 = match(j,nmatch(j))
@@ -440,5 +440,5 @@
        else
           ptop(1:npoints,ibox)=0.
-          do ilev=1,nlev
+          DO ilev=1,nlev
              where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0))
                 ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
@@ -458,6 +458,6 @@
     boxtau(1:npoints,1:ncol)  = output_missing_value
     boxptop(1:npoints,1:ncol) = output_missing_value
-    do ibox=1,ncol
-       do j=1,npoints 
+    DO ibox=1,ncol
+       DO j=1,npoints
           if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0.) then
              if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
@@ -532,10 +532,10 @@
     ! Reset frequencies
     !fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7)
-    do ilev=1,7
-       do ilev2=1,7
-          do j=1,npoints ! 
+    DO ilev=1,7
+       DO ilev2=1,7
+          DO j=1,npoints !
              if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then 
                 fq_isccp(j,ilev,ilev2)= 0.
-	     else 
+         else
                 fq_isccp(j,ilev,ilev2)= output_missing_value
              end if
@@ -559,5 +559,5 @@
     
     ! Compute column quantities and joint-histogram
-    do j=1,npoints 
+    DO j=1,npoints
        ! Subcolumns that are cloudy(true) and not(false)
        box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0.)
@@ -633,5 +633,5 @@
     
     varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
-    do j=1,dim2
+    DO j=1,dim2
        where(flag(:,j,:) .eq. 1)
           varOUT(:,j,:) = varIN2
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90	(revision 5158)
@@ -175,5 +175,5 @@
     ! PLANE PARRALLEL FIELDS
     ! ####################################################################################
-    do icol=1,ncolumns
+    DO icol=1,ncolumns
        ! #################################################################################
        ! *) Total Backscatter signal
@@ -200,5 +200,5 @@
     ! ####################################################################################
     if (lphaseoptics) then
-       do icol=1,ncolumns
+       DO icol=1,ncolumns
           ! #################################################################################
           ! *) Ice/Liq Perpendicular Backscatter signal
@@ -206,5 +206,5 @@
           ! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering 
           ! contribution (Cesana and Chepfer 2013, JGR)
-          do k=1,nlev
+          DO k=1,nlev
              ! Ice particles
              pnorm_perp_ice(1:npoints,icol,k) = Alpha * pnorm_ice(1:npoints,icol,k)
@@ -242,5 +242,5 @@
           
           ! Other layers
-          do k=2,nlev
+          DO k=2,nlev
              ! Optical thickness of layer k
              tautot_lay(1:npoints) = tautot(1:npoints,icol,k)-tautot(1:npoints,icol,k-1) 
@@ -398,5 +398,5 @@
     ! Compute LIDAR scattering ratio
     if (use_vgrid) then
-       do ic = 1, ncol
+       DO ic = 1, ncol
           pnorm_c = pnormFlip(:,ic,:)
           where ((pnorm_c .lt. xmax) .and. (betamolFlip(:,1,:) .lt. xmax) .and.          &
@@ -427,5 +427,5 @@
        endif
     else
-       do ic = 1, ncol
+       DO ic = 1, ncol
           pnorm_c = pnorm(:,ic,:)
           where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
@@ -458,6 +458,6 @@
     if (ok_lidar_cfad) then
        ! CFADs of subgrid-scale lidar scattering ratios
-       do i=1,Npoints
-          do j=1,llm
+       DO i=1,Npoints
+          DO j=1,llm
              cfad2(i,:,j) = hist1D(ncol,x3d(i,:,j),SR_BINS,histBsct)
           enddo
@@ -499,5 +499,5 @@
 
     ! Other layers
-    do k=2,nlev
+    DO k=2,nlev
        tautot_lay(:) = tau(:,k)-tau(:,k-1) 
        WHERE (tautot_lay(:) .gt. 0.)
@@ -527,5 +527,5 @@
     epsrealwp = epsilon(1._wp)
     beta(:,1) = pnorm(:,1) * (2._wp*tau(:,1))/(1._wp-exp(-2._wp*tau(:,1)))
-    do k=2,nlev
+    DO k=2,nlev
        tautot_lay(:) = tau(:,k)-tau(:,k-1)       
        WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. epsrealwp )
@@ -569,5 +569,5 @@
        zeta50    = -9.4776e-07_wp    ! 
        
-	! Inputs
+    ! Inputs
     integer,intent(in) :: &
        Npoints,  & ! Number of gridpoints
@@ -576,5 +576,5 @@
        Ncat,     & ! Number of cloud layer types
        Nphase      ! Number of cloud layer phase types
-	               ! [ice,liquid,undefined,false ice,false liquid,Percent of ice]
+                   ! [ice,liquid,undefined,false ice,false liquid,Percent of ice]
     real(wp),intent(in) :: &
        S_att,    & !
@@ -590,5 +590,5 @@
        pplay       ! Pressure
 
-	! Outputs
+    ! Outputs
     real(wp),intent(out),dimension(Npoints,Ntemp,5) :: &
        lidarcldtemp  ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
@@ -625,5 +625,5 @@
              
     ! ####################################################################################
-	! 1) Initialize    
+    ! 1) Initialize
     ! ####################################################################################
     lidarcld              = 0._wp
@@ -648,5 +648,5 @@
     ! 2) Cloud detection
     ! ####################################################################################
-    do k=1,Nlevels
+    DO k=1,Nlevels
        ! Cloud detection at subgrid-scale:
        where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
@@ -671,12 +671,12 @@
     cldlay   = 0._wp
     nsublay  = 0._wp
-    do k=1,Nlevels
-       do ic = 1, Ncolumns
-          do ip = 1, Npoints
+    DO k=1,Nlevels
+       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
+                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._wp
@@ -686,5 +686,5 @@
              
              if(cldy(ip,ic,k).eq.1.)then
-                do itemp=1,Ntemp 
+                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._wp
@@ -723,6 +723,6 @@
     cldlayer  = 0._wp
     nsublayer = 0._wp
-    do iz = 1, Ncat
-       do ic = 1, Ncolumns
+    DO iz = 1, Ncat
+       DO ic = 1, Ncolumns
           cldlayer(:,iz)  = cldlayer(:,iz)  + cldlay(:,ic,iz)
           nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
@@ -742,7 +742,7 @@
     ! 4.1) For Cloudy pixels with 8.16km < z < 19.2km
     ! ####################################################################################
-    do ncol=1,Ncolumns
-       do i=1,Npoints          
-          do nlev=1,23 ! from 19.2km until 8.16km
+    DO ncol=1,Ncolumns
+       DO i=1,Npoints
+          DO nlev=1,23 ! from 19.2km until 8.16km
                p1 = pplay(1,nlev)
 
@@ -856,5 +856,5 @@
           ! ##############################################################################
           toplvlsat = 0
-          do nlev=24,Nlevels! from 8.16km until 0km
+          DO nlev=24,Nlevels! from 8.16km until 0km
              p1 = pplay(i,nlev)
 
@@ -979,5 +979,5 @@
           ! ##############################################################################
           if(toplvlsat.ne.0) then
-             do nlev = toplvlsat,Nlevels
+             DO nlev = toplvlsat,Nlevels
                 p1 = pplay(i,nlev)
                 if(cldy(i,ncol,nlev).eq.1.)then
@@ -1031,7 +1031,7 @@
 
     ! Compute Phase low mid high cloud fractions
-    do iz = 1, Ncat
-       do i=1,Nphase-3
-          do ic = 1, Ncolumns
+    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)
@@ -1039,7 +1039,7 @@
        enddo
     enddo
-    do iz = 1, Ncat
-       do i=4,5
-          do ic = 1, Ncolumns
+    DO iz = 1, Ncat
+       DO i=4,5
+          DO ic = 1, Ncolumns
              cldlayerphase(:,iz,i) = cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
           enddo
@@ -1055,5 +1055,5 @@
     ENDWHERE
     
-    do i=1,Nphase-1
+    DO i=1,Nphase-1
        WHERE ( cldlayerphasesum(:,:).gt.0.0 )
           cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
@@ -1061,10 +1061,10 @@
     enddo
     
-    do i=1,Npoints
-       do iz=1,Ncat
+    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
+             DO ic=1,Nphase-3
                 checkcldlayerphase = checkcldlayerphase+cldlayerphase(i,iz,ic)
              enddo
@@ -1075,5 +1075,5 @@
     enddo
     
-    do i=1,Nphase-1
+    DO i=1,Nphase-1
        WHERE (nsublayer(:,:) .eq. 0.0)
           cldlayerphase(:,:,i) = undef
@@ -1082,8 +1082,8 @@
  
     ! Compute Phase 3D as a function of temperature
-    do nlev=1,Nlevels
-       do ncol=1,Ncolumns
-          do i=1,Npoints
-             do itemp=1,Ntemp
+    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
@@ -1105,6 +1105,6 @@
     
     ! Check temperature cloud fraction
-    do i=1,Npoints
-       do itemp=1,Ntemp
+    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
@@ -1124,5 +1124,5 @@
     ENDWHERE
     
-    do i=1,4
+    DO i=1,4
        WHERE(lidarcldtempind(:,:) .gt. 0.)
           lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
@@ -1191,5 +1191,5 @@
     ! 2) Cloud detection
     ! ####################################################################################
-    do k=1,Nlevels
+    DO k=1,Nlevels
        ! Cloud detection at subgrid-scale:
        where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
@@ -1210,7 +1210,7 @@
     ! 3) Grid-box 3D cloud fraction and layered cloud fractions(ISCCP pressure categories)
     ! ####################################################################################
-    do k=1,Nlevels
-       do ic = 1, Ncolumns
-          do ip = 1, Npoints
+    DO k=1,Nlevels
+       DO ic = 1, Ncolumns
+          DO ip = 1, Npoints
 
              iz=1
@@ -1244,6 +1244,6 @@
     cldlayer  = 0._wp
     nsublayer = 0._wp
-    do iz = 1, Ncat
-       do ic = 1, Ncolumns
+    DO iz = 1, Ncat
+       DO ic = 1, Ncolumns
           cldlayer(:,iz)  = cldlayer(:,iz)  + cldlay(:,ic,iz)
           nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
@@ -1272,5 +1272,5 @@
        eta = 0.6_wp            ! Multiple-scattering factor (Vaillant de Guelis et al. 2017a, AMT)
 
-	! Inputs
+    ! Inputs
     integer,intent(in) :: &
        Npoints,  & ! Number of gridpoints
@@ -1291,5 +1291,5 @@
        surfelev    ! Surface Elevation (SE)
 
-	! Outputs
+    ! Outputs
     real(wp),intent(out),dimension(Npoints,Nlevels,Ntype+1) :: &
        lidarcldtype   ! 3D OPAQ product fraction (opaque clouds, thin clouds, z_opaque, opacity)
@@ -1324,5 +1324,5 @@
 
     ! ####################################################################################
-	! 1) Initialize    
+    ! 1) Initialize
     ! ####################################################################################
     cldtype(:,:)          = 0._wp
@@ -1342,5 +1342,5 @@
     ! 2) Cloud detection and Fully attenuated layer detection
     ! ####################################################################################
-    do k=1,Nlevels
+    DO k=1,Nlevels
        ! Cloud detection at subgrid-scale:
        where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
@@ -1375,7 +1375,7 @@
     ! ####################################################################################
 
-    do k=1,Nlevels
-       do ic = 1, Ncolumns
-          do ip = 1, Npoints
+    DO k=1,Nlevels
+       DO ic = 1, Ncolumns
+          DO ip = 1, Npoints
 
              cldlay(ip,ic,1)   = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque cloud
@@ -1393,32 +1393,32 @@
 
 ! OPAQ variables
-     do ic = 1, Ncolumns
-        do ip = 1, Npoints
+     DO ic = 1, Ncolumns
+        DO ip = 1, Npoints
 
      ! Declaring non-opaque cloudy profiles as thin cloud profiles
-	   if ( cldlay(ip,ic,4).gt. 0. .and. cldlay(ip,ic,1) .eq. 0. ) then
-	      cldlay(ip,ic,2)  =  1._wp
- 	   endif
+       if ( cldlay(ip,ic,4).gt. 0. .and. cldlay(ip,ic,1) .eq. 0. ) then
+          cldlay(ip,ic,2)  =  1._wp
+        endif
 
      ! Filling in 3D and 2D variables
 
      ! Opaque cloud profiles
-	   if ( cldlay(ip,ic,1) .eq. 1. ) then
-	      zopac = 0._wp
-	      z_top = 0._wp
-	      do k=1,Nlevels-1
+       if ( cldlay(ip,ic,1) .eq. 1. ) then
+          zopac = 0._wp
+          z_top = 0._wp
+          DO k=1,Nlevels-1
      ! Declaring z_opaque altitude and opaque cloud fraction for 3D and 2D variables
      ! From SFC-2-TOA ( actually from vgrid_z(SFC+1) = vgrid_z(Nlevels-1) )
-	         if ( cldy(ip,ic,Nlevels-k) .eq. 1. .and. zopac .eq. 0. ) then
-		    lidarcldtype(ip,Nlevels-k + 1,3) = lidarcldtype(ip,Nlevels-k + 1,3) + 1._wp
-		    cldlay(ip,ic,3)                  = vgrid_z(Nlevels-k+1)      ! z_opaque altitude
-		    nsublay(ip,ic,3)                 = 1._wp
-		    zopac = Nlevels-k+1                        ! z_opaque vertical index on vgrid_z
-		 endif
-	         if ( cldy(ip,ic,Nlevels-k) .eq. 1. ) then
-		    lidarcldtype(ip,Nlevels-k ,1)    = lidarcldtype(ip,Nlevels-k ,1) + 1._wp
-		    z_top = Nlevels-k    ! top cloud layer vertical index on vgrid_z
+             if ( cldy(ip,ic,Nlevels-k) .eq. 1. .and. zopac .eq. 0. ) then
+            lidarcldtype(ip,Nlevels-k + 1,3) = lidarcldtype(ip,Nlevels-k + 1,3) + 1._wp
+            cldlay(ip,ic,3)                  = vgrid_z(Nlevels-k+1)      ! z_opaque altitude
+            nsublay(ip,ic,3)                 = 1._wp
+            zopac = Nlevels-k+1                        ! z_opaque vertical index on vgrid_z
+         endif
+             if ( cldy(ip,ic,Nlevels-k) .eq. 1. ) then
+            lidarcldtype(ip,Nlevels-k ,1)    = lidarcldtype(ip,Nlevels-k ,1) + 1._wp
+            z_top = Nlevels-k    ! top cloud layer vertical index on vgrid_z
                  endif
-	      enddo
+          enddo
      ! Summing opaque cloud mean temperatures and altitudes
      ! as defined in Vaillant de Guelis et al. 2017a, AMT
@@ -1432,44 +1432,44 @@
                  cldlay(ip,ic,1) = 0
               endif
-	   endif
+       endif
 
      ! Thin cloud profiles
-	   if ( cldlay(ip,ic,2) .eq. 1. ) then
-	      topcloud = 0._wp
-	      z_top = 0._wp
-	      z_base = 0._wp
-	      do k=1,Nlevels
+       if ( cldlay(ip,ic,2) .eq. 1. ) then
+          topcloud = 0._wp
+          z_top = 0._wp
+          z_base = 0._wp
+          DO k=1,Nlevels
      ! Declaring thin cloud fraction for 3D variable
      ! From TOA-2-SFC
                  if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 1. ) then
                     lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
-		    z_base = k ! bottom cloud layer
+            z_base = k ! bottom cloud layer
                  endif
-	         if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 0. ) then
+             if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 0. ) then
                     lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
-		    z_top = k  ! top cloud layer
-		    z_base = k ! bottom cloud layer
+            z_top = k  ! top cloud layer
+            z_base = k ! bottom cloud layer
                     topcloud = 1._wp
-		 endif
-	      enddo
+         endif
+          enddo
      ! Computing mean emissivity using layers below the bottom cloud layer to the surface
-     	      srmean = 0._wp
-	      srcount = 0._wp
-	      cloudemis = 0._wp
-     	      do k=z_base+1,Nlevels
-	         if (  (x(ip,ic,k) .gt. S_att_opaq) .and. (x(ip,ic,k) .lt. 1.0) .and. (x(ip,ic,k) .ne. undef)  ) then
-		    srmean = srmean + x(ip,ic,k)
-		    srcount = srcount + 1.
+               srmean = 0._wp
+          srcount = 0._wp
+          cloudemis = 0._wp
+               DO k=z_base+1,Nlevels
+             if (  (x(ip,ic,k) .gt. S_att_opaq) .and. (x(ip,ic,k) .lt. 1.0) .and. (x(ip,ic,k) .ne. undef)  ) then
+            srmean = srmean + x(ip,ic,k)
+            srcount = srcount + 1.
                  endif
-	      enddo
-	      ! If clear sky layers exist below bottom cloud layer
-	      if ( srcount .gt. 0. ) then
-	      	 trans2 = srmean/srcount              ! thin cloud transmittance**2
-	      	 tau_app = -(log(trans2))/2.          ! apparent cloud optical depth
-	      	 tau_vis = tau_app/eta                ! cloud visible optical depth (multiple scat.)
-	      	 tau_ir = tau_vis/2.                  ! approx. relation between visible and IR ODs
-	      	 cloudemis = 1. - exp(-tau_ir)        ! no diffusion in IR considered : emis = 1-T
-		 count_emis(ip) = count_emis(ip) + 1.
-	      endif
+          enddo
+          ! If clear sky layers exist below bottom cloud layer
+          if ( srcount .gt. 0. ) then
+               trans2 = srmean/srcount              ! thin cloud transmittance**2
+               tau_app = -(log(trans2))/2.          ! apparent cloud optical depth
+               tau_vis = tau_app/eta                ! cloud visible optical depth (multiple scat.)
+               tau_ir = tau_vis/2.                  ! approx. relation between visible and IR ODs
+               cloudemis = 1. - exp(-tau_ir)        ! no diffusion in IR considered : emis = 1-T
+         count_emis(ip) = count_emis(ip) + 1.
+          endif
      ! Summing thin cloud mean temperatures and altitudes
      ! as defined in Vaillant de Guelis et al. 2017a, AMT
@@ -1500,18 +1500,18 @@
     ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=1) to SFC(k=Nlevels)
        lidarcldtype(:,1,4) = lidarcldtype(:,1,3) !top layer equal to 3D z_opaque fraction
-    do ip = 1, Npoints
-     	do k = 2, Nlevels
+    DO ip = 1, Npoints
+         DO k = 2, Nlevels
             if ( (lidarcldtype(ip,k,3) .ne. undef) .and. (lidarcldtype(ip,k-1,4) .ne. undef) ) then
-	        lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4)
-	    else
-	        lidarcldtype(ip,k,4) = undef
-	    endif
-	enddo
+            lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4)
+        else
+            lidarcldtype(ip,k,4) = undef
+        endif
     enddo
+    enddo
 
     ! Layered cloud types (opaque, thin and z_opaque 2D variables)
 
-    do iz = 1, Ntype
-       do ic = 1, Ncolumns
+    DO iz = 1, Ntype
+       DO ic = 1, Ncolumns
           cldtype(:,iz)  = cldtype(:,iz)  + cldlay(:,ic,iz)
           nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90	(revision 5158)
@@ -15,7 +15,7 @@
                               numMODISReffIceBins,reffICE_binCenters, &
                               numMODISReffLiqBins, reffLIQ_binCenters, &
-			      numISCCPTauBins,numISCCPPresBins, &
-			      numMISRTauBins,numMISRHgtBins, &
-			      numModisTauBins,numMODISPresBins
+                  numISCCPTauBins,numISCCPPresBins, &
+                  numMISRTauBins,numMISRHgtBins, &
+                  numModisTauBins,numMODISPresBins
 
   implicit none
@@ -297,7 +297,7 @@
     type(cosp_optical_inputs),intent(inout) :: y
 
-    if (allocated(y%tau_067))       	  deallocate(y%tau_067)
-    if (allocated(y%emiss_11))      	  deallocate(y%emiss_11)
-    if (allocated(y%frac_out))    	  deallocate(y%frac_out)
+    if (allocated(y%tau_067))             deallocate(y%tau_067)
+    if (allocated(y%emiss_11))            deallocate(y%emiss_11)
+    if (allocated(y%frac_out))          deallocate(y%frac_out)
     if (allocated(y%beta_mol_calipso))    deallocate(y%beta_mol_calipso)
     if (allocated(y%tau_mol_calipso))     deallocate(y%tau_mol_calipso)
@@ -308,12 +308,12 @@
     if (allocated(y%tautot_ice_calipso))  deallocate(y%tautot_ice_calipso)
     if (allocated(y%tautot_liq_calipso))  deallocate(y%tautot_liq_calipso)
-    if (allocated(y%tautot_S_liq))	  deallocate(y%tautot_S_liq)
-    if (allocated(y%tautot_S_ice))    	  deallocate(y%tautot_S_ice)
-    if (allocated(y%z_vol_cloudsat))  	  deallocate(y%z_vol_cloudsat)
-    if (allocated(y%kr_vol_cloudsat)) 	  deallocate(y%kr_vol_cloudsat)
-    if (allocated(y%g_vol_cloudsat)) 	  deallocate(y%g_vol_cloudsat)
-    if (allocated(y%asym))            	  deallocate(y%asym)
-    if (allocated(y%ss_alb))         	  deallocate(y%ss_alb)
-    if (allocated(y%fracLiq))       	  deallocate(y%fracLiq)
+    if (allocated(y%tautot_S_liq))      deallocate(y%tautot_S_liq)
+    if (allocated(y%tautot_S_ice))          deallocate(y%tautot_S_ice)
+    if (allocated(y%z_vol_cloudsat))        deallocate(y%z_vol_cloudsat)
+    if (allocated(y%kr_vol_cloudsat))       deallocate(y%kr_vol_cloudsat)
+    if (allocated(y%g_vol_cloudsat))       deallocate(y%g_vol_cloudsat)
+    if (allocated(y%asym))                  deallocate(y%asym)
+    if (allocated(y%ss_alb))               deallocate(y%ss_alb)
+    if (allocated(y%fracLiq))             deallocate(y%fracLiq)
     if (allocated(y%beta_mol_grLidar532)) deallocate(y%beta_mol_grLidar532)
     if (allocated(y%betatot_grLidar532))  deallocate(y%betatot_grLidar532)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90	(revision 5158)
@@ -53,6 +53,6 @@
 ! mr_ozone,                             !Concentration ozone (Kg/Kg)
 ! dem_s                                 !Cloud optical emissivity
-! dtau_s               			!Cloud optical thickness
-! emsfc_lw = 1.        			!Surface emissivity dans radlwsw.F90
+! dtau_s                           !Cloud optical thickness
+! emsfc_lw = 1.                    !Surface emissivity dans radlwsw.F90
 
 
@@ -95,40 +95,40 @@
 
   ! Local variables
-  character(len=64),PARAMETER	:: cosp_input_nl  = 'cospv2_input_nl.txt'
-  character(len=64),PARAMETER  	:: cosp_output_nl = 'cospv2_output_nl.txt'
-
-  integer, save 	       	:: isccp_topheight, isccp_topheight_direction, overlap
-  integer, save  	       	:: Ncolumns    	    	! Number of subcolumns in SCOPS
-  integer, save 	       	:: Npoints            	! Number of gridpoints
+  character(len=64),PARAMETER    :: cosp_input_nl  = 'cospv2_input_nl.txt'
+  character(len=64),PARAMETER      :: cosp_output_nl = 'cospv2_output_nl.txt'
+
+  integer, save                :: isccp_topheight, isccp_topheight_direction, overlap
+  integer, save                 :: Ncolumns                ! Number of subcolumns in SCOPS
+  integer, save                :: Npoints                ! Number of gridpoints
 !$OMP THREADPRIVATE(Npoints)
-  integer, save			:: Nlevels		! Number of model vertical levels
-  integer  			:: Nptslmdz, Nlevlmdz  	! Nb de points issus de physiq.F
-  integer, save 		:: Npoints_it   	! Max number of gridpoints to be 
-  	   			   			! processed in one iteration
-  type(cosp_config), save 	:: cfg          	! Variable qui contient les cles
-  		     		   			! logiques des simulateurs et des
-							! diagnostics, definie dans:
-							! lmdz_cosp_construct_destroy_mod
+  integer, save            :: Nlevels        ! Number of model vertical levels
+  integer              :: Nptslmdz, Nlevlmdz      ! Nb de points issus de physiq.F
+  integer, save         :: Npoints_it       ! Max number of gridpoints to be
+                                    ! processed in one iteration
+  type(cosp_config), save     :: cfg              ! Variable qui contient les cles
+                                      ! logiques des simulateurs et des
+                            ! diagnostics, definie dans:
+                            ! lmdz_cosp_construct_destroy_mod
 !$OMP THREADPRIVATE(cfg)
 
-  integer			:: t0, t1, count_rate, count_max
-  real(wp), save 		:: cloudsat_radar_freq, cloudsat_k2, rttov_ZenAng, co2, &
-  	    			   ch4, n2o, co, emsfc_lw
+  integer            :: t0, t1, count_rate, count_max
+  real(wp), save         :: cloudsat_radar_freq, cloudsat_k2, rttov_ZenAng, co2, &
+                         ch4, n2o, co, emsfc_lw
 !$OMP THREADPRIVATE(emsfc_lw)
 
-  integer, dimension(RTTOV_MAX_CHANNELS), save	 :: rttov_Channels
+  integer, dimension(RTTOV_MAX_CHANNELS), save     :: rttov_Channels
   real(wp), dimension(RTTOV_MAX_CHANNELS), save  :: rttov_Surfem
-  integer, save 			   	 :: surface_radar, use_mie_tables, &
-  	   					    cloudsat_use_gas_abs, cloudsat_do_ray, &
-						    melt_lay
-  integer, save 				 :: lidar_ice_type
-  integer, save 				 :: rttov_platform, rttov_satellite, &
-  	   					    rttov_Instrument, rttov_Nchannels
-  logical, save 				 :: use_vgrid_in, csat_vgrid_in, &
-  	   					    use_precipitation_fluxes
+  integer, save                     :: surface_radar, use_mie_tables, &
+                                 cloudsat_use_gas_abs, cloudsat_do_ray, &
+                            melt_lay
+  integer, save                  :: lidar_ice_type
+  integer, save                  :: rttov_platform, rttov_satellite, &
+                                 rttov_Instrument, rttov_Nchannels
+  logical, save                  :: use_vgrid_in, csat_vgrid_in, &
+                                 use_precipitation_fluxes
 
 ! Declaration necessaires pour les sorties IOIPSL
-  real        	:: ecrit_day, ecrit_hf, ecrit_mth, missing_val
-  logical 	:: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml
+  REAL            :: ecrit_day, ecrit_hf, ecrit_mth, missing_val
+  logical     :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml
   logical, save :: debut_cosp=.true.
 !$OMP THREADPRIVATE(debut_cosp)
@@ -142,19 +142,19 @@
 
 !-----------------------------  Input variables from LMDZ-GCM  -------------------------------
-  integer                             	     :: overlaplmdz   ! overlap type: 1=max, 
-  					       		      ! 2=rand, 3=max/rand
-  real, dimension(Nptslmdz,Nlevlmdz) 	     :: phi, p, ph, T, sh, rh, tca, cca, mr_lsliq,   &
-  					        mr_lsice, mr_ccliq, mr_ccice, fl_lsrain,     & 
-                                       	        fl_lssnow, fl_ccrain, fl_ccsnow, fl_lsgrpl,  &
-                                       	        zlev, zlev_half, mr_ozone, radliq, radice,   &
-					        dtau_s, dem_s, dtau_c, dem_c, ref_liq, ref_ice
-  real, dimension(Nptslmdz,Nlevlmdz) 	     :: fl_lsrainI, fl_lssnowI, fl_ccrainI, fl_ccsnowI
-  real, dimension(Nptslmdz)        	     :: lon, lat, skt, fracTerLic, u_wind, v_wind, &
-  					        phis, sunlit
-  real, dimension(Nptslmdz)        	     :: land ! variables intermediaire pour masque TerLic
-  real, dimension(Nlevlmdz)        	     :: presnivs
-  integer                         	     :: itap, k, ip
-  real                            	     :: dtime, freq_cosp
-  real, dimension(2)               	     :: time_bnds
+  integer                                      :: overlaplmdz   ! overlap type: 1=max,
+                                           ! 2=rand, 3=max/rand
+  real, dimension(Nptslmdz,Nlevlmdz)          :: phi, p, ph, T, sh, rh, tca, cca, mr_lsliq,   &
+                              mr_lsice, mr_ccliq, mr_ccice, fl_lsrain,     &
+                                                   fl_lssnow, fl_ccrain, fl_ccsnow, fl_lsgrpl,  &
+                                                   zlev, zlev_half, mr_ozone, radliq, radice,   &
+                            dtau_s, dem_s, dtau_c, dem_c, ref_liq, ref_ice
+  real, dimension(Nptslmdz,Nlevlmdz)          :: fl_lsrainI, fl_lssnowI, fl_ccrainI, fl_ccsnowI
+  real, dimension(Nptslmdz)                 :: lon, lat, skt, fracTerLic, u_wind, v_wind, &
+                              phis, sunlit
+  real, dimension(Nptslmdz)                 :: land ! variables intermediaire pour masque TerLic
+  real, dimension(Nlevlmdz)                 :: presnivs
+  integer                                  :: itap, k, ip
+  REAL                                     :: dtime, freq_cosp
+  real, dimension(2)                        :: time_bnds
 
   double precision                           :: d_dtime
@@ -169,20 +169,20 @@
   logical :: &
        Lsingle     = .true.,  & ! True if using MMF_v3_single_moment CLOUDSAT
-       		     	      	! microphysical scheme (default)
+                                  ! microphysical scheme (default)
        Ldouble     = .false.    ! True if using MMF_v3.5_two_moment CLOUDSAT 
-       		     		! microphysical scheme
-  type(size_distribution), save   	     :: sd            ! Hydrometeor description
+                            ! microphysical scheme
+  type(size_distribution), save            :: sd            ! Hydrometeor description
 !$OMP THREADPRIVATE(sd)
-  type(radar_cfg), save			     :: rcfg_cloudsat ! Radar configuration
+  type(radar_cfg), save                 :: rcfg_cloudsat ! Radar configuration
 !$OMP THREADPRIVATE(rcfg_cloudsat)
-  real, dimension(Nptslmdz,Nlevlmdz,N_HYDRO) :: Reff	      ! Liquid and Ice particles
-  					     		      ! effective radius
-  type(cosp_outputs)			     :: cospOUT       ! COSP simulator outputs
-  type(cosp_optical_inputs) 		     :: cospIN        ! COSP optical (or derived?)
-  					     		      ! fields needed by simulators
-  type(cosp_column_inputs) 		     :: cospstateIN   ! COSP model fields needed
-  					     		      ! by simulators
-  character(len=256), dimension(100)	     :: cosp_status
-  character(len=64), save      		     :: cloudsat_micro_scheme
+  real, dimension(Nptslmdz,Nlevlmdz,N_HYDRO) :: Reff          ! Liquid and Ice particles
+                                         ! effective radius
+  type(cosp_outputs)                 :: cospOUT       ! COSP simulator outputs
+  type(cosp_optical_inputs)              :: cospIN        ! COSP optical (or derived?)
+                                         ! fields needed by simulators
+  type(cosp_column_inputs)              :: cospstateIN   ! COSP model fields needed
+                                         ! by simulators
+  character(len=256), dimension(100)         :: cosp_status
+  character(len=64), save                   :: cloudsat_micro_scheme
 
   ! Indices to address arrays of LS and CONV hydrometeors
@@ -271,5 +271,5 @@
   print*,' Cles des differents simulateurs cosp a itap :',itap
   print*,'cfg%Lcloudsat, cfg%Lcalipso, cfg%LgrLidar532, cfg%Latlid, cfg%Lparasol, &
-	cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov', &
+    cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov', &
         cfg%Lcloudsat, cfg%Lcalipso, cfg%LgrLidar532, cfg%Latlid, cfg%Lparasol, &
         cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov
@@ -340,7 +340,7 @@
 print*,' Cles des differents simulateurs cosp a itap :',itap
 print*,'cfg%Lcloudsat, cfg%Lcalipso, cfg%LgrLidar532, cfg%Latlid, cfg%Lparasol, &
-	cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov', &
+    cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov', &
         cfg%Lcloudsat, cfg%Lcalipso, cfg%LgrLidar532, cfg%Latlid, cfg%Lparasol, &
-	cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov
+    cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov
 
   endif !(itap.gt.1).and.(first_write)
@@ -369,6 +369,6 @@
 
         zlev_half(:,1) = phis(:)/9.81
-        do k = 2, Nlevels
-          do ip = 1, Npoints
+        DO k = 2, Nlevels
+          DO ip = 1, Npoints
            zlev_half(ip,k) = phi(ip,k)/9.81 + &
                (phi(ip,k)-phi(ip,k-1))/9.81 * (ph(ip,k)-p(ip,k)) / (p(ip,k)-p(ip,k-1))
@@ -383,5 +383,5 @@
 
 ! 3) Masque terre/mer a partir de la variable fracTerLic
-        do ip = 1, Npoints
+        DO ip = 1, Npoints
           if (fracTerLic(ip).ge.0.5) then
              land(ip) = 1.
@@ -445,19 +445,19 @@
     call construct_cospstateIN(Npoints,Nlevels,0,cospstateIN)
 
-    cospstateIN%lat                 		        = lat(1:Npoints)
-    cospstateIN%lon                     		= lon(1:Npoints) 
-    cospstateIN%at              	      		= t(1:Npoints,Nlevels:1:-1) 
-    cospstateIN%qv              	     		= sh(1:Npoints,Nlevels:1:-1)
-    cospstateIN%o3              	      		= mr_ozone(1:Npoints,Nlevels:1:-1)  
-    cospstateIN%sunlit          	      		= sunlit(1:Npoints)
-    cospstateIN%skt             	       		= skt(1:Npoints)
-    cospstateIN%land            	     		= land(1:Npoints)
-    cospstateIN%surfelev            	     		= zlev_half(1:Npoints,1)
-    cospstateIN%pfull           	      		= p(1:Npoints,Nlevels:1:-1)
-    cospstateIN%phalf(1:Npoints,1)	                = 0._wp
-    cospstateIN%phalf(1:Npoints,2:Nlevels+1)	        = ph(1:Npoints,Nlevels:1:-1)  
-    cospstateIN%hgt_matrix                  		= zlev(1:Npoints,Nlevels:1:-1) 
-    cospstateIN%hgt_matrix_half(1:Npoints,Nlevels+1) 	= 0._wp
-    cospstateIN%hgt_matrix_half(1:Npoints,1:Nlevels) 	= zlev_half(1:Npoints,Nlevels:1:-1) 
+    cospstateIN%lat                                 = lat(1:Npoints)
+    cospstateIN%lon                             = lon(1:Npoints)
+    cospstateIN%at                                = t(1:Npoints,Nlevels:1:-1)
+    cospstateIN%qv                               = sh(1:Npoints,Nlevels:1:-1)
+    cospstateIN%o3                                = mr_ozone(1:Npoints,Nlevels:1:-1)
+    cospstateIN%sunlit                            = sunlit(1:Npoints)
+    cospstateIN%skt                                = skt(1:Npoints)
+    cospstateIN%land                             = land(1:Npoints)
+    cospstateIN%surfelev                             = zlev_half(1:Npoints,1)
+    cospstateIN%pfull                             = p(1:Npoints,Nlevels:1:-1)
+    cospstateIN%phalf(1:Npoints,1)                    = 0._wp
+    cospstateIN%phalf(1:Npoints,2:Nlevels+1)            = ph(1:Npoints,Nlevels:1:-1)
+    cospstateIN%hgt_matrix                          = zlev(1:Npoints,Nlevels:1:-1)
+    cospstateIN%hgt_matrix_half(1:Npoints,Nlevels+1)     = 0._wp
+    cospstateIN%hgt_matrix_half(1:Npoints,1:Nlevels)     = zlev_half(1:Npoints,Nlevels:1:-1)
 
 
@@ -483,21 +483,21 @@
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-    call subsample_and_optics(cfg, Npoints, Nlevels, Ncolumns, N_HYDRO,overlap,	&
-        		      use_precipitation_fluxes, lidar_ice_type, sd, &
-			      tca(1:Npoints,Nlevels:1:-1), cca(1:Npoints,Nlevels:1:-1), &
-           		      fl_lsrain(1:Npoints,Nlevels:1:-1), &
-			      fl_lssnow(1:Npoints,Nlevels:1:-1), &
-			      fl_lsgrpl(1:Npoints,Nlevels:1:-1), &
-			      fl_ccrain(1:Npoints,Nlevels:1:-1), &
-          		      fl_ccsnow(1:Npoints,Nlevels:1:-1), &
-			      mr_lsliq(1:Npoints,Nlevels:1:-1), &
-			      mr_lsice(1:Npoints,Nlevels:1:-1), &
-			      mr_ccliq(1:Npoints,Nlevels:1:-1), &
-			      mr_ccice(1:Npoints,Nlevels:1:-1), &
-             		      Reff(1:Npoints,Nlevels:1:-1,:), &
-			      dtau_c(1:Npoints,Nlevels:1:-1), &
-			      dtau_s(1:Npoints,Nlevels:1:-1), &
-			      dem_c(1:Npoints,Nlevels:1:-1), &
-			      dem_s(1:Npoints,Nlevels:1:-1), cospstateIN, cospIN)
+    call subsample_and_optics(cfg, Npoints, Nlevels, Ncolumns, N_HYDRO,overlap,    &
+                      use_precipitation_fluxes, lidar_ice_type, sd, &
+                  tca(1:Npoints,Nlevels:1:-1), cca(1:Npoints,Nlevels:1:-1), &
+                         fl_lsrain(1:Npoints,Nlevels:1:-1), &
+                  fl_lssnow(1:Npoints,Nlevels:1:-1), &
+                  fl_lsgrpl(1:Npoints,Nlevels:1:-1), &
+                  fl_ccrain(1:Npoints,Nlevels:1:-1), &
+                        fl_ccsnow(1:Npoints,Nlevels:1:-1), &
+                  mr_lsliq(1:Npoints,Nlevels:1:-1), &
+                  mr_lsice(1:Npoints,Nlevels:1:-1), &
+                  mr_ccliq(1:Npoints,Nlevels:1:-1), &
+                  mr_ccice(1:Npoints,Nlevels:1:-1), &
+                           Reff(1:Npoints,Nlevels:1:-1,:), &
+                  dtau_c(1:Npoints,Nlevels:1:-1), &
+                  dtau_s(1:Npoints,Nlevels:1:-1), &
+                  dem_c(1:Npoints,Nlevels:1:-1), &
+                  dem_s(1:Npoints,Nlevels:1:-1), cospstateIN, cospIN)
 
 
@@ -531,5 +531,5 @@
                                    ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml, &
                                    ecrit_mth, ecrit_day, ecrit_hf, use_vgrid_in,         &
-			           niv_sorties, vgrid_z_in, zlev(1,:))
+                       niv_sorties, vgrid_z_in, zlev(1,:))
 
       !$OMP END MASTER
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_mod.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_mod.F90	(revision 5158)
@@ -88,5 +88,5 @@
          "clcalipsoliq", "CALIPSO Liq-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_clcalipsoun = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
-         "clcalipsoun", "CALIPSO Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))	 
+         "clcalipsoun", "CALIPSO Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_clcalipsotmpice = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
          "clcalipsotmpice", "CALIPSO Ice-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
@@ -109,5 +109,5 @@
          "clcalipsothin", "CALIPSO Thin profile Cloud Fraction", "%", (/ ('', i=1, 3) /))      
   TYPE(ctrl_outcosp), SAVE :: o_clcalipsozopaque = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
-         "clcalipsozopaque", "CALIPSO z_opaque Fraction", "%", (/ ('', i=1, 3) /))	   
+         "clcalipsozopaque", "CALIPSO z_opaque Fraction", "%", (/ ('', i=1, 3) /))
   TYPE(ctrl_outcosp), SAVE :: o_clcalipsoopacity = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
          "clcalipsoopacity", "CALIPSO opacity Fraction", "%", (/ ('', i=1, 3) /))
@@ -311,5 +311,5 @@
                               ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml,  &
                               ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, Nlvgrid, vgrid_z_loc, &
-			      vgrid_mz)
+                  vgrid_mz)
 
   use mod_cosp_config, only : CLOUDSAT_DBZE_BINS, SR_BINS, CLOUDSAT_CFAD_ZE_MIN, PARASOL_NREFL, &
@@ -319,5 +319,5 @@
                               numMODISReffIceBins,reffICE_binCenters, &
                               numMODISReffLiqBins, reffLIQ_binCenters, pres_binCenters, &
-			      cloudsat_binCenters, calipso_binCenters
+                  cloudsat_binCenters, calipso_binCenters
 
   USE iophy
@@ -335,5 +335,5 @@
   real,dimension(Nlevlmdz) :: presnivs, vgrid_mz
   real,dimension(Nlvgrid)  :: vgrid_z_loc
-  real                     :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth 
+  REAL                     :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth
   logical                  :: use_vgrid
   logical                  :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml
@@ -341,5 +341,5 @@
 !!! Variables locales
   integer                   :: idayref, iff, ii
-  real                      :: zjulian,zjulian_start
+  REAL                      :: zjulian,zjulian_start
   real(wp),dimension(Ncolumns)  :: column_ax
   CHARACTER(LEN=20), DIMENSION(3)  :: chfreq = (/ '1day', '1d  ', '3h  ' /)            
@@ -359,5 +359,5 @@
 
 !! Definition valeurs axes
-    do ii=1,Ncolumns
+    DO ii=1,Ncolumns
       column_ax(ii) = real(ii)
     enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90	(revision 5158)
@@ -46,5 +46,5 @@
 !!! Variables d'entree
   integer               :: itap, Nlevlmdz, Ncolumns, Npoints, Nlvgrid
-  real                  :: freq_COSP, dtime, missing_val, missing_cosp
+  REAL                  :: freq_COSP, dtime, missing_val, missing_cosp
   type(cosp_config)     :: cfg     ! Control outputs
   type(cosp_outputs)    :: &
@@ -257,7 +257,7 @@
      where(cospOUT%calipso_cfad_sr == R_UNDEF) cospOUT%calipso_cfad_sr = missing_val
     
-     do icl=1,SR_BINS
-      do k=1,Nlvgrid
-       do ip=1,Npoints
+     DO icl=1,SR_BINS
+      DO k=1,Nlvgrid
+       DO ip=1,Npoints
         tmp_fi4da_cfadL(ip,k,icl)=cospOUT%calipso_cfad_sr(ip,icl,k)
        enddo
@@ -312,7 +312,7 @@
    if (cfg%LcfadLidarsr532gr) then
      where(cospOUT%grLidar532_cfad_sr == R_UNDEF) cospOUT%grLidar532_cfad_sr = missing_val
-     do icl=1,SR_BINS
-      do k=1,Nlvgrid
-       do ip=1,Npoints
+     DO icl=1,SR_BINS
+      DO k=1,Nlvgrid
+       DO ip=1,Npoints
          tmp_fi4da_cfadLgr(ip,k,icl)=cospOUT%grLidar532_cfad_sr(ip,icl,k)
        enddo
@@ -360,7 +360,7 @@
    if (cfg%LcfadLidarsr355) then
      where(cospOUT%atlid_cfad_sr == R_UNDEF) cospOUT%atlid_cfad_sr = missing_val
-     do icl=1,SR_BINS
-      do k=1,Nlvgrid
-       do ip=1,Npoints
+     DO icl=1,SR_BINS
+      DO k=1,Nlvgrid
+       DO ip=1,Npoints
           tmp_fi4da_cfadLatlid(ip,k,icl)=cospOUT%atlid_cfad_sr(ip,icl,k)
        enddo
@@ -456,7 +456,7 @@
    if (cfg%LcfadDbze94)  then
      where(cospOUT%cloudsat_cfad_ze == R_UNDEF) cospOUT%cloudsat_cfad_ze = missing_val  
-     do icl=1,CLOUDSAT_DBZE_BINS
-      do k=1,Nlvgrid
-       do ip=1,Npoints
+     DO icl=1,CLOUDSAT_DBZE_BINS
+      DO k=1,Nlvgrid
+       DO ip=1,Npoints
          tmp_fi4da_cfadR(ip,k,icl)=cospOUT%cloudsat_cfad_ze(ip,icl,k)
        enddo
@@ -531,7 +531,7 @@
    where(cospOUT%misr_cldarea == R_UNDEF) cospOUT%misr_cldarea = missing_val
 
-   do icl=1,numMISRHgtBins
-      do k=1,Nlvgrid
-       do ip=1,Npoints    
+   DO icl=1,numMISRHgtBins
+      DO k=1,Nlvgrid
+       DO ip=1,Npoints
       tmp_fi4da_misr(ip,icl,k)=cospOUT%misr_fq(ip,k,icl)
        enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90	(revision 5158)
@@ -29,5 +29,5 @@
          Lcloudsat,        & ! CLOUDSAT simulator on/off switch 
          Lcalipso,         & ! CALIPSO simulator on/off switch 
-	 LgrLidar532,      & ! GROUND LIDAR simulator on/off switch
+     LgrLidar532,      & ! GROUND LIDAR simulator on/off switch
          Latlid,           & ! ATLID simulator on/off switch
          Lisccp,           & ! ISCCP simulator on/off switch
@@ -161,5 +161,5 @@
   integer :: i
                 
-   do i=1,107
+   DO i=1,107
       cfg%out_list(i)=''
    enddo
@@ -295,5 +295,5 @@
   integer :: i
 
-   do i=1,107
+   DO i=1,107
       cfg%out_list(i)=''
    enddo
@@ -490,5 +490,5 @@
                        Lptradarflag8,Lptradarflag9,Lradarpia
 
-  do i=1,107
+  DO i=1,107
     cfg%out_list(i)=''
   enddo
@@ -1062,5 +1062,5 @@
  IF (using_xios) THEN
 
-  do i=1,107
+  DO i=1,107
     cfg%out_list(i)=''
   enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90	(revision 5158)
@@ -142,7 +142,7 @@
        frac_cv(1:nPoints,1:nLevels) = 0._wp
        prec_cv(1:nPoints,1:nLevels) = 0._wp
-       do j=1,nPoints
-          do k=1,nLevels
-             do i=1,nColumns
+       DO j=1,nPoints
+          DO k=1,nLevels
+             DO i=1,nColumns
                 if (cospIN%frac_out(j,i,k)  .eq. 1)  frac_ls(j,k) = frac_ls(j,k)+1._wp
                 if (cospIN%frac_out(j,i,k)  .eq. 2)  frac_cv(j,k) = frac_cv(j,k)+1._wp
@@ -171,5 +171,5 @@
        Reff(:,:,:,:)     = 0._wp
        Np(:,:,:,:)       = 0._wp
-       do k=1,nColumns
+       DO k=1,nColumns
           ! Subcolumn cloud fraction
           column_frac_out = cospIN%frac_out(:,k,:)
@@ -214,6 +214,6 @@
        fl_ccrain(:,:) = 0._wp
        fl_ccsnow(:,:) = 0._wp
-       do k=1,nLevels
-          do j=1,nPoints
+       DO k=1,nLevels
+          DO j=1,nPoints
              ! In-cloud mixing ratios.
              if (frac_ls(j,k) .ne. 0.) then
@@ -359,6 +359,6 @@
        allocate(g_vol(nPoints,nLevels))
        g_vol(:,:)=0._wp
-       do i=1,nPoints
-          do j=1,nLevels
+       DO i=1,nPoints
+          DO j=1,nLevels
              if (cospIN%rcfg_cloudsat%use_gas_abs == 1 .or. (cospIN%rcfg_cloudsat%use_gas_abs == 2 .and. j .eq. 1)) then
                 g_vol(i,j) = gases(cospstateIN%pfull(i,j), cospstateIN%at(i,j),cospstateIN%qv(i,j),cospIN%rcfg_cloudsat%freq)
@@ -371,5 +371,5 @@
        allocate(fracPrecipIce(nPoints,nColumns,nLevels))
        fracPrecipIce(:,:,:) = 0._wp
-       do k=1,nColumns
+       DO k=1,nColumns
           call quickbeam_optics(sd, cospIN%rcfg_cloudsat, nPoints, nLevels, R_UNDEF,  &
                mr_hydro(:,k,:,1:nHydro)*1000._wp, Reff(:,k,:,1:nHydro)*1.e6_wp,&
@@ -406,6 +406,6 @@
 
        ! For near-surface diagnostics, we only need the frozen fraction at one layer.
-       do i=1,nPoints                                                                       !PREC_BUG
-	  cospIN%fracPrecipIce(i,:) = fracPrecipIce_statGrid(i,:,cloudsat_preclvl_index(i)) !PREC_BUG
+       DO i=1,nPoints                                                                       !PREC_BUG
+      cospIN%fracPrecipIce(i,:) = fracPrecipIce_statGrid(i,:,cloudsat_preclvl_index(i)) !PREC_BUG
        enddo                                                                                !PREC_BUG
 !       cospIN%fracPrecipIce(:,:) = fracPrecipIce_statGrid(:,:,cloudsat_preclvl)      !!! ORIGINAL
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/math_lib.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/math_lib.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/math_lib.F90	(revision 5158)
@@ -96,5 +96,5 @@
     else
        sumo = 0._wp
-       do j=i1,i2
+       DO j=i1,i2
           deltah = abs(s(i1+1)-s(i1))
           sumo = sumo + f(j)*deltah
@@ -204,5 +204,5 @@
     end if
     
-    do i = 2, ntab
+    DO i = 2, ntab
        if ( xtab(i) <= xtab(i-1) ) then
           lerror = .true.
@@ -239,5 +239,5 @@
     ihi = ntab
     
-    do i = 1, ntab
+    DO i = 1, ntab
        if ( a <= xtab(i) ) then
           exit
@@ -249,5 +249,5 @@
     ilo = min ( ilo, ntab - 1 )
     
-    do i = 1, ntab
+    DO i = 1, ntab
        if ( xtab(i) <= b ) then
           exit
@@ -263,5 +263,5 @@
 !ds    sum1 = 0.0D+00
     
-    do i = ilo, ihi
+    DO i = ilo, ihi
        
        x1 = xtab(i-1)
@@ -371,5 +371,5 @@
           ga=1._wp
           m1=x-1
-          do k=2,m1
+          DO k=2,m1
              ga=ga*k
           enddo
@@ -382,5 +382,5 @@
           m=int(z)
           r=1._wp
-          do k=1,m
+          DO k=1,m
              r=r*(z-k)
           enddo
@@ -390,5 +390,5 @@
        endif
        gr=g(26)
-       do k=25,1,-1
+       DO k=25,1,-1
           gr=gr*z+g(k)
        enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/modis_simulator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/modis_simulator.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/modis_simulator.F90	(revision 5158)
@@ -163,5 +163,5 @@
     cloudMask = retrievedTau(1:nSubCols) >= min_OpticalThickness
     
-    do i = 1, nSubCols
+    DO i = 1, nSubCols
        if(cloudMask(i)) then 
           ! ##################################################################################
@@ -380,5 +380,5 @@
     reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
     reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
-    do j=1,nPoints
+    DO j=1,nPoints
 
        ! Fill clear and optically thin subcolumns with fill
@@ -439,5 +439,5 @@
     !   layers and use the trapezoidal rule.
     totalTau = 0._wp; totalProduct = 0._wp
-    do i = 2, size(tauIncrement)
+    DO i = 2, size(tauIncrement)
       if(totalTau + tauIncrement(i) > tauLimit) then 
         deltaX = tauLimit - totalTau
@@ -479,5 +479,5 @@
     ! Find the extinction-weighted value of f(tau), assuming constant f within each layer
     totalTau = 0._wp; totalProduct = 0._wp
-    do i = 1, size(tauIncrement)
+    DO i = 1, size(tauIncrement)
       if(totalTau + tauIncrement(i) > tauLimit) then 
         deltaX       = tauLimit - totalTau
@@ -712,5 +712,5 @@
     cloudMask(1:nLevels) = tau(1:nLevels) > 0. 
     cloudIndicies = pack((/ (i, i = 1, nLevels) /), mask = cloudMask) 
-    do i = 1, size(cloudIndicies)
+    DO i = 1, size(cloudIndicies)
        call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
     end do
@@ -892,5 +892,5 @@
     Tran_cumulative(1) = Tran(1)    
     
-    do i=2, npts
+    DO i=2, npts
        ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
        Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1._wp - Refl_cumulative(i-1) * Refl(i))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90	(revision 5158)
@@ -539,5 +539,5 @@
     if (alam < cutice) then
        ! Region from 0.045 microns to 167.0 microns - no temperature depend
-       do i=2,nwl
+       DO i=2,nwl
           if(alam < wl(i)) continue
        enddo
@@ -557,10 +557,10 @@
        if(tk > temref(1)) tk=temref(1)
        if(tk < temref(4)) tk=temref(4)
-       do i=2,4
+       DO i=2,4
           if(tk.ge.temref(i)) go to 12
        enddo
 12     lt1 = i
        lt2 = i-1
-       do i=2,nwlt
+       DO i=2,nwlt
           if(alam.le.wlt(i)) go to 14
        enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90	(revision 5158)
@@ -159,11 +159,11 @@
     ! Compute grid-box averaged Parasol reflectances
     parasolrefl(:,:) = 0._wp
-    do k = 1, nrefl
-       do ic = 1, ncol
+    DO k = 1, nrefl
+       DO ic = 1, ncol
           parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
        enddo
     enddo
     
-    do k = 1, nrefl
+    DO k = 1, nrefl
        parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
        ! if land=1 -> parasolrefl=R_UNDEF
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90	(revision 5158)
@@ -66,7 +66,7 @@
       if (cv_col .eq. 0) cv_col=1
  
-      do ilev=1,nlev
-        do ibox=1,ncol
-          do j=1,npoints 
+      DO ilev=1,nlev
+        DO ibox=1,ncol
+          DO j=1,npoints
             prec_frac(j,ibox,ilev) = 0
           enddo
@@ -74,11 +74,11 @@
       enddo
       
-      do j=1,npoints
-       do ibox=1,ncol
+      DO j=1,npoints
+       DO ibox=1,ncol
         frac_out_ls(j,ibox)=0
         frac_out_cv(j,ibox)=0
         flag_ls=0
         flag_cv=0
-        do ilev=1,nlev
+        DO ilev=1,nlev
           if (frac_out(j,ibox,ilev) .eq. 1) then 
             flag_ls=1
@@ -98,10 +98,10 @@
 
 !      initialize the top layer      
-       do j=1,npoints
+       DO j=1,npoints
         flag_ls=0
         flag_cv=0
     
         if (ls_p_rate(j,1) .gt. 0.) then 
-            do ibox=1,ncol ! possibility ONE
+            DO ibox=1,ncol ! possibility ONE
                 if (frac_out(j,ibox,1) .eq. 1) then 
                     prec_frac(j,ibox,1) = 1
@@ -110,5 +110,5 @@
             enddo ! loop over ncol
             if (flag_ls .eq. 0) then ! possibility THREE
-                do ibox=1,ncol
+                DO ibox=1,ncol
                     if (frac_out(j,ibox,2) .eq. 1) then 
                         prec_frac(j,ibox,1) = 1
@@ -118,5 +118,5 @@
             endif
         if (flag_ls .eq. 0) then ! possibility Four
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out_ls(j,ibox) .eq. 1) then 
             prec_frac(j,ibox,1) = 1
@@ -126,5 +126,5 @@
         endif
         if (flag_ls .eq. 0) then ! possibility Five
-        do ibox=1,ncol
+        DO ibox=1,ncol
     !     prec_frac(j,1:ncol,1) = 1
         prec_frac(j,ibox,1) = 1
@@ -135,5 +135,5 @@
      
         if (cv_p_rate(j,1) .gt. 0.) then 
-         do ibox=1,ncol ! possibility ONE
+         DO ibox=1,ncol ! possibility ONE
           if (frac_out(j,ibox,1) .eq. 2) then 
            if (prec_frac(j,ibox,1) .eq. 0) then
@@ -146,5 +146,5 @@
         enddo ! loop over ncol
         if (flag_cv .eq. 0) then ! possibility THREE
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out(j,ibox,2) .eq. 2) then 
                 if (prec_frac(j,ibox,1) .eq. 0) then
@@ -158,5 +158,5 @@
         endif
         if (flag_cv .eq. 0) then ! possibility Four
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out_cv(j,ibox) .eq. 1) then 
                 if (prec_frac(j,ibox,1) .eq. 0) then
@@ -170,5 +170,5 @@
         endif
         if (flag_cv .eq. 0) then  ! possibility Five
-        do ibox=1,cv_col
+        DO ibox=1,cv_col
                 if (prec_frac(j,ibox,1) .eq. 0) then
             prec_frac(j,ibox,1) = 2
@@ -187,11 +187,11 @@
 
 !     working on the levels from top to surface
-      do ilev=2,nlev
-       do j=1,npoints
+      DO ilev=2,nlev
+       DO j=1,npoints
         flag_ls=0
         flag_cv=0
     
         if (ls_p_rate(j,ilev) .gt. 0.) then 
-         do ibox=1,ncol ! possibility ONE&TWO
+         DO ibox=1,ncol ! possibility ONE&TWO
           if ((frac_out(j,ibox,ilev) .eq. 1) .or. ((prec_frac(j,ibox,ilev-1) .eq. 1)     &
             .or. (prec_frac(j,ibox,ilev-1) .eq. 3))) then 
@@ -201,5 +201,5 @@
         enddo ! loop over ncol
         if ((flag_ls .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out(j,ibox,ilev+1) .eq. 1) then 
             prec_frac(j,ibox,ilev) = 1
@@ -209,5 +209,5 @@
         endif
         if (flag_ls .eq. 0) then ! possibility Four
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out_ls(j,ibox) .eq. 1) then 
             prec_frac(j,ibox,ilev) = 1
@@ -217,5 +217,5 @@
         endif
         if (flag_ls .eq. 0) then ! possibility Five
-        do ibox=1,ncol
+        DO ibox=1,ncol
 !     prec_frac(j,1:ncol,ilev) = 1
         prec_frac(j,ibox,ilev) = 1
@@ -225,5 +225,5 @@
     
         if (cv_p_rate(j,ilev) .gt. 0.) then 
-         do ibox=1,ncol ! possibility ONE&TWO
+         DO ibox=1,ncol ! possibility ONE&TWO
           if ((frac_out(j,ibox,ilev) .eq. 2) .or. ((prec_frac(j,ibox,ilev-1) .eq. 2)     &
             .or. (prec_frac(j,ibox,ilev-1) .eq. 3))) then 
@@ -237,5 +237,5 @@
        enddo ! loop over ncol
         if ((flag_cv .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out(j,ibox,ilev+1) .eq. 2) then 
                 if (prec_frac(j,ibox,ilev) .eq. 0) then
@@ -249,5 +249,5 @@
         endif
         if (flag_cv .eq. 0) then ! possibility Four
-        do ibox=1,ncol
+        DO ibox=1,ncol
         if (frac_out_cv(j,ibox) .eq. 1) then 
                 if (prec_frac(j,ibox,ilev) .eq. 0) then
@@ -261,5 +261,5 @@
         endif
         if (flag_cv .eq. 0) then  ! possibility Five 
-        do ibox=1,cv_col
+        DO ibox=1,cv_col
                 if (prec_frac(j,ibox,ilev) .eq. 0) then
             prec_frac(j,ibox,ilev) = 2
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90	(revision 5158)
@@ -150,7 +150,7 @@
        d_gate     = -1
     endif
-    do k=start_gate,end_gate,d_gate
+    DO k=start_gate,end_gate,d_gate
        ! Loop over each profile (nprof)
-       do pr=1,nprof
+       DO pr=1,nprof
           ! Attenuation due to hydrometeors between radar and volume
           
@@ -284,6 +284,6 @@
           
           ! Effective reflectivity histogram
-          do i=1,Npoints
-             do j=1,llm
+          DO i=1,Npoints
+             DO j=1,llm
                 cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_toti(i,:,j),DBZE_BINS,cloudsat_histRef)
              enddo
@@ -301,6 +301,6 @@
        else
           ! Effective reflectivity histogram
-          do i=1,Npoints
-             do j=1,llm
+          DO i=1,Npoints
+             DO j=1,llm
                 cfad_ze(i,:,j) = hist1D(Ncolumns,Ze_tot(i,:,j),DBZE_BINS,cloudsat_histRef)
              enddo
@@ -395,10 +395,10 @@
     ! SUBCOLUMN processing
     ! ######################################################################################
-    do i=1, Npoints
+    DO i=1, Npoints
 
        cloudsat_preclvl = cloudsat_preclvl_index(i) !PREC_BUG
 !       print*, 'i, surfelev(i), cloudsat_preclvl : ', i, surfelev(i), cloudsat_preclvl !PREC_BUG
 
-       do pr=1,Ncolumns
+       DO pr=1,Ncolumns
           ! 1) Compute the PIA in all profiles containing hydrometeors
           if ( (Ze_non_out(i,pr,cloudsat_preclvl).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then
@@ -523,7 +523,7 @@
 
     ! Aggregate subcolumns
-    do i=1,Npoints
+    DO i=1,Npoints
        ! Gridmean precipitation fraction for each precipitation type
-       do k=1,nCloudsatPrecipClass
+       DO k=1,nCloudsatPrecipClass
           if (any(cloudsat_pflag(i,:) .eq. k-1)) then
              cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq. k-1)
@@ -571,7 +571,7 @@
             trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
        
-       do i=1,maxhclass
-          do j=1,mt_ntt
-             do k=1,nRe_types
+       DO i=1,maxhclass
+          DO j=1,mt_ntt
+             DO k=1,nRe_types
                 ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
                 read(12,rec=ind) rcfg%Z_scale_flag(i,j,k), &
@@ -607,7 +607,7 @@
          trim(rcfg%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
     
-    do i=1,maxhclass
-       do j=1,mt_ntt
-          do k=1,nRe_types
+    DO i=1,maxhclass
+       DO j=1,mt_ntt
+          DO k=1,nRe_types
              ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
              if(.not.LUT_file_exists .or. rcfg%Z_scale_added_flag(i,j,k)) then
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90	(revision 5158)
@@ -71,5 +71,5 @@
     mt_ttl = (/ ((j-1)*5-60 + 273.15, j = 1, cnt_liq) /)
     D(1) = dmin
-    do j=2,nd
+    DO j=2,nd
        D(j) = D(j-1)*exp((log(dmax)-log(dmin))/(nd-1))
     enddo
@@ -140,9 +140,9 @@
     kr_vol   = 0._wp
 
-    do k=1,ngate       ! Loop over each profile (nprof)
-       do pr=1,nprof
+    DO k=1,ngate       ! Loop over each profile (nprof)
+       DO pr=1,nprof
           ! Determine if hydrometeor(s) present in volume
           hydro = .false.
-          do j=1,rcfg%nhclass
+          DO j=1,rcfg%nhclass
              if ((hm_matrix(pr,k,j) > 1E-12) .and. (sd%dtype(j) > 0)) then
                 hydro = .true.
@@ -157,5 +157,5 @@
              
              ! Loop over hydrometeor type
-             do tp=1,rcfg%nhclass
+             DO tp=1,rcfg%nhclass
                 Re_internal = re_matrix(pr,k,tp)
 
@@ -342,5 +342,5 @@
     
     where(kr_vol(:,:) <= EPSILON(kr_vol)) 
-       ! Volume is hydrometeor-free	
+       ! Volume is hydrometeor-free
        !z_vol(:,:)  = undef
        z_ray(:,:)  = undef
@@ -794,5 +794,5 @@
        lidx = infind(D,dmin)
        uidx = infind(D,dmax)    
-       do k=lidx,uidx
+       DO k=lidx,uidx
           N(k) = (ahp*(D(k)*1E-3)**bhp) * 1E-12    
        enddo
@@ -980,5 +980,5 @@
        sizep = (pi*D0)/wl
        dqv(1) = 0._wp
-       do i=1,nsizes
+       DO i=1,nsizes
           call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), &
                dg, xs1, xs2, dph, err)
@@ -1183,5 +1183,5 @@
     sumo = 0._wp
     aux1 = 1.1_wp*e_th
-    do i=1,nbands_o2
+    DO i=1,nbands_o2
        aux2   = f/v0(i)
        aux3   = v0(i)-f
@@ -1207,5 +1207,5 @@
     sumo = 0._wp
     aux1 = 4.8_wp*e_th
-    do i=1,nbands_h2o
+    DO i=1,nbands_h2o
        aux2    = f/v1(i)
        aux3    = v1(i)-f
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90	(revision 5157)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90	(revision 5158)
@@ -94,5 +94,5 @@
     if (ncolprint.ne.0) then
        write (6,'(a)') 'frac_out_pp_rev:'
-       do j=1,npoints,1000
+       DO j=1,npoints,1000
           write(6,'(a10)') 'j='
           write(6,'(8I10)') j
@@ -104,5 +104,5 @@
     if (ncolprint.ne.0) then
        write (6,'(a)') 'last_frac_pp:'
-       do j=1,npoints,1000
+       DO j=1,npoints,1000
           write(6,'(a10)') 'j='
           write(6,'(8I10)') j
@@ -139,5 +139,5 @@
           IF (ncolprint.ne.0) then
              write (6,'(a)') 'threshold_nsf2:'
-             do j=1,npoints,1000
+             DO j=1,npoints,1000
                 write(6,'(a10)') 'j='
                 write(6,'(8I10)') j
@@ -156,5 +156,5 @@
           !maxocc(1:npoints,ibox) = merge(1,0,boxpos(1:npoints,ibox) .le. conv(1:npoints,ilev))
           !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox))
-          do j=1,npoints
+          DO j=1,npoints
              if (boxpos(j,ibox).le.conv(j,ilev)) then
                 maxocc(j,ibox) = 1
@@ -214,5 +214,5 @@
        ! Set last_frac to tca at this level, so as to be tca from last level next time round
        if (ncolprint.ne.0) then
-          do j=1,npoints ,1000
+          DO j=1,npoints ,1000
              write(6,'(a10)') 'j='
              write(6,'(8I10)') j
