Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/MISR_simulator.F90	(revision 5082)
@@ -84,5 +84,5 @@
        do ilev=1,nlev
           ! Define location of "layer top"
-          if(ilev.eq.1 .or. ilev.eq.nlev) then
+          if(ilev==1 .or. ilev==nlev) then
              ztest=zfull(j,ilev)
           else
@@ -94,5 +94,5 @@
           iMISR_ztop=2
           do loop=2,numMISRHgtBins
-             if ( ztest .gt. 1000*misr_histHgt(loop+1) ) then
+             if ( ztest > 1000*misr_histHgt(loop+1) ) then
                 iMISR_ztop=loop+1
              endif
@@ -110,5 +110,5 @@
           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
+             if(thres_crossed_MISR == 0 .and. dtau(j,ibox,ilev) > 0.) then
                 ! First encountered a "cloud"
                 thres_crossed_MISR = 1  
@@ -116,6 +116,6 @@
              endif
 
-             if( thres_crossed_MISR .lt. 99 .and. thres_crossed_MISR .gt. 0 ) then
-                if( dtau(j,ibox,ilev) .eq. 0.) then
+             if( thres_crossed_MISR < 99 .and. thres_crossed_MISR > 0 ) then
+                if( dtau(j,ibox,ilev) == 0.) then
                    ! We have come to the end of the current cloud layer without yet 
                    ! selecting a CTH boundary. Restart cloud tau counter 
@@ -129,6 +129,6 @@
                 ! current layer cloud top to the current level then MISR will like 
                 ! see a top below the top of the current layer.
-                if( dtau(j,ibox,ilev).gt.0 .and. (cloud_dtau-dtau(j,ibox,ilev)) .lt. 1) then
-                   if(dtau(j,ibox,ilev) .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then
+                if( dtau(j,ibox,ilev)>0 .and. (cloud_dtau-dtau(j,ibox,ilev)) < 1) then
+                   if(dtau(j,ibox,ilev) < 1 .or. ilev==1 .or. ilev==nlev) then
                       ! MISR will likely penetrate to some point within this layer ... the middle
                       MISR_penetration_height=zfull(j,ilev)
@@ -142,5 +142,5 @@
                 
                 ! Check for a distinctive water layer
-                if(dtau(j,ibox,ilev) .gt. 1 .and. at(j,ilev) .gt. 273 ) then
+                if(dtau(j,ibox,ilev) > 1 .and. at(j,ilev) > 273 ) then
                    ! Must be a water cloud, take this as CTH level
                    thres_crossed_MISR=99
@@ -149,5 +149,5 @@
                 ! If the total column optical depth is "large" than MISR can't see
                 ! anything else. Set current point as CTH level
-                if(sum(dtau(j,ibox,1:ilev)) .gt. 5) then
+                if(sum(dtau(j,ibox,1:ilev)) > 5) then
                    thres_crossed_MISR=99           
                 endif
@@ -157,12 +157,12 @@
           ! Check to see if there was a cloud for which we didn't 
           ! set a MISR cloud top boundary
-          if( thres_crossed_MISR .eq. 1) then
+          if( thres_crossed_MISR == 1) then
              ! If the cloud has a total optical depth of greater
              ! than ~ 0.5 MISR will still likely pick up this cloud
              ! with a height near the true cloud top
              ! otherwise there should be no CTH
-             if(sum(dtau(j,ibox,1:nlev)) .gt. 0.5) then
+             if(sum(dtau(j,ibox,1:nlev)) > 0.5) then
                 ! keep MISR detected CTH
-             elseif(sum(dtau(j,ibox,1:nlev)) .gt. 0.2) then
+             elseif(sum(dtau(j,ibox,1:nlev)) > 0.2) then
                 ! MISR may detect but wont likley have a good height
                 box_MISR_ztop(j,ibox)=-1
@@ -215,5 +215,5 @@
     ! Fill dark scenes 
     do j=1,numMISRHgtBins
-       where(sunlit .ne. 1) dist_model_layertops(1:npoints,j) = R_UNDEF
+       where(sunlit /= 1) dist_model_layertops(1:npoints,j) = R_UNDEF
     enddo
 
@@ -257,12 +257,12 @@
 
        ! Subcolumns that are cloudy(true) and not(false)
-       box_cloudy(1:ncol) = merge(.true.,.false.,tau(j,1:ncol) .gt. tauchk)
+       box_cloudy(1:ncol) = merge(.true.,.false.,tau(j,1:ncol) > tauchk)
 
        ! Fill optically thin clouds with fill value
        where(.not. box_cloudy(1:ncol)) tauWRK(j,1:ncol)  = -999._wp
-       where(box_MISR_ztopWRK(j,1:ncol) .eq. 0) box_MISR_ztopWRK(j,1:ncol)=-999._wp
+       where(box_MISR_ztopWRK(j,1:ncol) == 0) box_MISR_ztopWRK(j,1:ncol)=-999._wp
 
        ! Compute joint histogram and column quantities for points that are sunlit and cloudy
-       if (sunlit(j) .eq. 1) then 
+       if (sunlit(j) == 1) then
           ! Joint histogram
           call hist2D(tauWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol),ncol,misr_histTau,numMISRTauBins,&
@@ -272,10 +272,10 @@
 
           ! Column cloud area
-          MISR_cldarea(j)=real(count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.))/ncol
+          MISR_cldarea(j)=real(count(box_MISR_ztopWRK(j,1:ncol) /= -999.))/ncol
 
           ! Column cloud-top height
-          if ( count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.) .ne. 0 ) then
-             MISR_mean_ztop(j) = sum(box_MISR_ztopWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol) .ne. -999.)/ &
-                  count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.)
+          if ( count(box_MISR_ztopWRK(j,1:ncol) /= -999.) /= 0 ) then
+             MISR_mean_ztop(j) = sum(box_MISR_ztopWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol) /= -999.)/ &
+                  count(box_MISR_ztopWRK(j,1:ncol) /= -999.)
           else
              MISR_mean_ztop(j) = R_UNDEF
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp.F90	(revision 5082)
@@ -480,5 +480,5 @@
 
     ! Set flag to deallocate rttov types (only done on final call to simulator)
-    if (size(cospOUT%isccp_meantb) .eq. stop_idx) lrttov_cleanUp = .true.
+    if (size(cospOUT%isccp_meantb) == stop_idx) lrttov_cleanUp = .true.
 
     ! ISCCP column
@@ -687,10 +687,10 @@
        modisIN%w0        => cospIN%ss_alb
        modisIN%Nsunlit   = count(cospgridIN%sunlit > 0)
-       if (modisIN%Nsunlit .gt. 0) then
+       if (modisIN%Nsunlit > 0) then
           allocate(modisIN%sunlit(modisIN%Nsunlit),modisIN%pres(modisIN%Nsunlit,cospIN%Nlevels+1))
           modisIN%sunlit    = pack((/ (i, i = 1, Npoints ) /),mask = cospgridIN%sunlit > 0)
           modisIN%pres      = cospgridIN%phalf(int(modisIN%sunlit(:)),:)
        endif
-       if (count(cospgridIN%sunlit <= 0) .gt. 0) then
+       if (count(cospgridIN%sunlit <= 0) > 0) then
           allocate(modisIN%notSunlit(count(cospgridIN%sunlit <= 0)))
           modisIN%notSunlit = pack((/ (i, i = 1, Npoints ) /),mask = .not. cospgridIN%sunlit > 0)
@@ -2430,5 +2430,5 @@
     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if (any([Lisccp_subcolumn, Lisccp_column, Lmisr_subcolumn, Lmisr_column, Lmodis_subcolumn, Lmodis_column])) then
-       if (any(cospgridIN%sunlit .lt. 0)) then
+       if (any(cospgridIN%sunlit < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%sunlit contains values out of range (0 or 1)'
@@ -2498,5 +2498,5 @@
          Lcalipso_column, Lcloudsat_column, Lradar_lidar_tcc,Llidar_only_freq_cloud, &
          Lcloudsat_tcc, Lcloudsat_tcc2])) then
-       if (any(cospgridIN%at .lt. 0)) then
+       if (any(cospgridIN%at < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%at contains values out of range (at<0), expected units (K)'
@@ -2546,5 +2546,5 @@
     endif
     if (any([Lisccp_subcolumn, Lisccp_column, Lrttov_column])) then
-       if (any(cospgridIN%pfull .lt. 0)) then
+       if (any(cospgridIN%pfull < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%pfull contains values out of range'
@@ -2566,5 +2566,5 @@
     if (any([Lisccp_subcolumn,Lisccp_column,Lmodis_subcolumn,Lmodis_column,Lcalipso_column,Lrttov_column,&
              LgrLidar532_column,Latlid_column])) then
-       if (any(cospgridIN%phalf .lt. 0)) then
+       if (any(cospgridIN%phalf < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%phalf contains values out of range'
@@ -2648,5 +2648,5 @@
     endif
     if (any([Lisccp_subcolumn,Lisccp_column,Lrttov_column])) then
-       if (any(cospgridIN%qv .lt. 0)) then
+       if (any(cospgridIN%qv < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%qv contains values out of range'
@@ -2668,5 +2668,5 @@
     if (any([Lmisr_subcolumn,Lmisr_column,Lcloudsat_subcolumn,Lcloudsat_column,Lcalipso_column,Lradar_lidar_tcc,&
          Llidar_only_freq_cloud,LgrLidar532_column,Latlid_column,Lcloudsat_tcc, Lcloudsat_tcc2])) then
-       if (any(cospgridIN%hgt_matrix .lt. -300)) then
+       if (any(cospgridIN%hgt_matrix < -300)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%hgt_matrix contains values out of range'
@@ -2714,5 +2714,5 @@
     if (any([Lrttov_column,Lcloudsat_column,Lcalipso_column,Lradar_lidar_tcc,Llidar_only_freq_cloud, &
              LgrLidar532_column, Latlid_column, Lcloudsat_tcc, Lcloudsat_tcc2])) then
-       if (any(cospgridIN%hgt_matrix_half .lt. -300)) then
+       if (any(cospgridIN%hgt_matrix_half < -300)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%hgt_matrix_half contains values out of range'
@@ -2753,5 +2753,5 @@
     endif
     if (any([Lrttov_column,Lcalipso_column,Lparasol_column])) then
-       if (any(cospgridIN%land .lt. 0)) then
+       if (any(cospgridIN%land < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%land contains values out of range'
@@ -2776,5 +2776,5 @@
     endif
     if (any([Lisccp_subcolumn,Lisccp_column,Lrttov_column])) then
-       if (any(cospgridIN%skt .lt. 0)) then
+       if (any(cospgridIN%skt < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%skt contains values out of range'
@@ -2797,5 +2797,5 @@
     ! RTTOV Inputs
     if (Lrttov_column) then
-       if (cospgridIN%zenang .lt. -90. .OR. cospgridIN%zenang .gt. 90) then
+       if (cospgridIN%zenang < -90. .OR. cospgridIN%zenang > 90) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%zenang contains values out of range'
@@ -2803,5 +2803,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (cospgridIN%co2 .lt. 0) then
+       if (cospgridIN%co2 < 0) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%co2 contains values out of range'
@@ -2809,5 +2809,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (cospgridIN%ch4 .lt. 0) then
+       if (cospgridIN%ch4 < 0) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%ch4 contains values out of range'
@@ -2815,5 +2815,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (cospgridIN%n2o .lt. 0) then
+       if (cospgridIN%n2o < 0) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%n2o contains values out of range'
@@ -2821,5 +2821,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (cospgridIN%co.lt. 0) then
+       if (cospgridIN%co< 0) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%co contains values out of range'
@@ -2827,5 +2827,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (any(cospgridIN%o3 .lt. 0)) then
+       if (any(cospgridIN%o3 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%o3 contains values out of range'
@@ -2833,5 +2833,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (any(cospgridIN%emis_sfc .lt. 0. .OR. cospgridIN%emis_sfc .gt. 1)) then
+       if (any(cospgridIN%emis_sfc < 0. .OR. cospgridIN%emis_sfc > 1)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospgridIN%emis_sfc contains values out of range'
@@ -2839,5 +2839,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (any(cospgridIN%u_sfc .lt. -100. .OR. cospgridIN%u_sfc .gt. 100.)) then
+       if (any(cospgridIN%u_sfc < -100. .OR. cospgridIN%u_sfc > 100.)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%u_sfc contains values out of range'
@@ -2845,5 +2845,5 @@
           Lrttov_column = .false.
        endif
-       if (any(cospgridIN%v_sfc .lt. -100. .OR. cospgridIN%v_sfc .gt. 100.)) then
+       if (any(cospgridIN%v_sfc < -100. .OR. cospgridIN%v_sfc > 100.)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%v_sfc contains values out of range'
@@ -2851,5 +2851,5 @@
           if (associated(cospOUT%rttov_tbs)) cospOUT%rttov_tbs(:,:) = R_UNDEF
        endif
-       if (any(cospgridIN%lat .lt. -90 .OR. cospgridIN%lat .gt. 90)) then
+       if (any(cospgridIN%lat < -90 .OR. cospgridIN%lat > 90)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%lat contains values out of range'
@@ -2861,5 +2861,5 @@
     ! COSP_INPUTS
     if (any([Lisccp_subcolumn,Lisccp_column])) then
-       if (cospIN%emsfc_lw .lt. 0. .OR. cospIN%emsfc_lw .gt. 1.) then
+       if (cospIN%emsfc_lw < 0. .OR. cospIN%emsfc_lw > 1.) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%emsfc_lw contains values out of range'
@@ -2878,5 +2878,5 @@
     endif
     if (any([Lisccp_subcolumn,Lisccp_column,Lmisr_subcolumn,Lmisr_column,Lmodis_subcolumn,Lmodis_column])) then
-       if (any(cospIN%tau_067 .lt. 0)) then
+       if (any(cospIN%tau_067 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tau_067 contains values out of range'
@@ -2943,5 +2943,5 @@
     endif
     if (any([Lisccp_subcolumn,Lisccp_column])) then
-       if (any(cospIN%emiss_11 .lt. 0. .OR. cospIN%emiss_11 .gt. 1)) then
+       if (any(cospIN%emiss_11 < 0. .OR. cospIN%emiss_11 > 1)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%emiss_11 contains values out of range'
@@ -2960,5 +2960,5 @@
     endif
     if (any([Lmodis_subcolumn,Lmodis_column])) then
-       if (any(cospIN%asym .lt. -1. .OR. cospIN%asym .gt. 1)) then
+       if (any(cospIN%asym < -1. .OR. cospIN%asym > 1)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%asym contains values out of range'
@@ -3006,5 +3006,5 @@
                cospOUT%modis_Optical_Thickness_vs_ReffLIQ(:,:,:)            = R_UNDEF
        endif
-       if (any(cospIN%ss_alb .lt. 0 .OR. cospIN%ss_alb .gt. 1)) then
+       if (any(cospIN%ss_alb < 0 .OR. cospIN%ss_alb > 1)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%ss_alb contains values out of range'
@@ -3054,5 +3054,5 @@
     endif
     if (any([Latlid_subcolumn,Latlid_column])) then
-       if (any(cospIN%betatot_atlid .lt. 0)) then
+       if (any(cospIN%betatot_atlid < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%betatot_atlid contains values out of range'
@@ -3065,5 +3065,5 @@
           if (associated(cospOUT%atlid_beta_mol))      cospOUT%atlid_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%beta_mol_atlid .lt. 0)) then
+       if (any(cospIN%beta_mol_atlid < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%beta_mol_atlid contains values out of range'
@@ -3076,5 +3076,5 @@
           if (associated(cospOUT%atlid_beta_mol))      cospOUT%atlid_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%tautot_atlid .lt. 0)) then
+       if (any(cospIN%tautot_atlid < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_atlid contains values out of range'
@@ -3087,5 +3087,5 @@
           if (associated(cospOUT%atlid_beta_mol))      cospOUT%atlid_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%tau_mol_atlid .lt. 0)) then
+       if (any(cospIN%tau_mol_atlid < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tau_mol_atlid contains values out of range'
@@ -3101,5 +3101,5 @@
     
     if (any([LgrLidar532_subcolumn,LgrLidar532_column])) then
-       if (any(cospIN%betatot_grLidar532 .lt. 0)) then
+       if (any(cospIN%betatot_grLidar532 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%betatot_grLidar532 contains values out of range'
@@ -3112,5 +3112,5 @@
           if (associated(cospOUT%grLidar532_beta_mol))      cospOUT%grLidar532_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%beta_mol_grLidar532 .lt. 0)) then
+       if (any(cospIN%beta_mol_grLidar532 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%beta_mol_grLidar532 contains values out of range'
@@ -3123,5 +3123,5 @@
           if (associated(cospOUT%grLidar532_beta_mol))      cospOUT%grLidar532_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%tautot_grLidar532 .lt. 0)) then
+       if (any(cospIN%tautot_grLidar532 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_grLidar532 contains values out of range'
@@ -3134,5 +3134,5 @@
           if (associated(cospOUT%grLidar532_beta_mol))      cospOUT%grLidar532_beta_mol(:,:)   = R_UNDEF
        endif
-       if (any(cospIN%tau_mol_grLidar532 .lt. 0)) then
+       if (any(cospIN%tau_mol_grLidar532 < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tau_mol_grLidar532 contains values out of range'
@@ -3148,5 +3148,5 @@
 
     if (any([Lcalipso_subcolumn,Lcalipso_column])) then
-       if (any(cospIN%betatot_calipso .lt. 0)) then
+       if (any(cospIN%betatot_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%betatot_calipso contains values out of range'
@@ -3167,5 +3167,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%betatot_liq_calipso .lt. 0)) then
+       if (any(cospIN%betatot_liq_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = ('ERROR: COSP input variable: cospIN%betatot_liq_calipso contains values out of range')
@@ -3186,5 +3186,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%betatot_ice_calipso .lt. 0)) then
+       if (any(cospIN%betatot_ice_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%betatot_ice_calipso contains values out of range'
@@ -3205,5 +3205,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%tautot_calipso .lt. 0)) then
+       if (any(cospIN%tautot_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_calipso contains values out of range'
@@ -3224,5 +3224,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%tautot_liq_calipso .lt. 0)) then
+       if (any(cospIN%tautot_liq_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = ('ERROR: COSP input variable: cospIN%tautot_liq_calipso contains values out of range')
@@ -3243,5 +3243,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%tautot_ice_calipso .lt. 0)) then
+       if (any(cospIN%tautot_ice_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_ice_calipso contains values out of range'
@@ -3262,5 +3262,5 @@
           if (associated(cospOUT%calipso_cldthinemis))   cospOUT%calipso_cldthinemis(:)       = R_UNDEF
        endif
-       if (any(cospIN%tau_mol_calipso .lt. 0)) then
+       if (any(cospIN%tau_mol_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tau_mol_calipso contains values out of range'
@@ -3284,5 +3284,5 @@
     if (any([Lcalipso_subcolumn,Lcalipso_column,Lcloudsat_column,Lradar_lidar_tcc,       &
         Llidar_only_freq_cloud, Lcloudsat_tcc, Lcloudsat_tcc2])) then
-       if (any(cospIN%beta_mol_calipso .lt. 0)) then
+       if (any(cospIN%beta_mol_calipso < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%beta_mol_calipso contains values out of range'
@@ -3315,5 +3315,5 @@
     endif
     if (any([Lparasol_subcolumn,Lparasol_column])) then
-       if (any(cospIN%tautot_S_liq .lt. 0)) then
+       if (any(cospIN%tautot_S_liq < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_S_liq contains values out of range'
@@ -3323,5 +3323,5 @@
           if (associated(cospOUT%parasolGrid_refl)) cospOUT%parasolGrid_refl(:,:)  = R_UNDEF
        endif
-       if (any(cospIN%tautot_S_ice .lt. 0)) then
+       if (any(cospIN%tautot_S_ice < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%tautot_S_ice contains values out of range'
@@ -3334,5 +3334,5 @@
     if (any([Lcloudsat_subcolumn,Lcloudsat_column,Lradar_lidar_tcc,Llidar_only_freq_cloud, &
         Lcloudsat_tcc, Lcloudsat_tcc2])) then
-       if (any(cospIN%z_vol_cloudsat .lt. 0)) then
+       if (any(cospIN%z_vol_cloudsat < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%z_vol_cloudsat contains values out of range'
@@ -3350,5 +3350,5 @@
           if (associated(cospOUT%cloudsat_tcc2)) cospOUT%cloudsat_tcc2(:) = R_UNDEF
        endif
-       if (any(cospIN%kr_vol_cloudsat .lt. 0)) then
+       if (any(cospIN%kr_vol_cloudsat < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%kr_vol_cloudsat contains values out of range'
@@ -3366,5 +3366,5 @@
           if (associated(cospOUT%cloudsat_tcc2)) cospOUT%cloudsat_tcc2(:) = R_UNDEF
        endif
-       if (any(cospIN%g_vol_cloudsat .lt. 0)) then
+       if (any(cospIN%g_vol_cloudsat < 0)) then
           nError=nError+1
           errorMessage(nError) = 'ERROR: COSP input variable: cospIN%g_vol_cloudsat contains values out of range'
@@ -3389,13 +3389,13 @@
     ! ISCCP
     if (Lisccp_subcolumn .or. Lisccp_column) then
-       if (size(cospIN%frac_out,1)  .ne. cospIN%Npoints .OR. &
-           size(cospIN%tau_067,1)   .ne. cospIN%Npoints .OR. &
-           size(cospIN%emiss_11,1)  .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%skt)     .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%qv,1)    .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%at,1)    .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%phalf,1) .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%sunlit)  .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%pfull,1) .ne. cospIN%Npoints) then
+       if (size(cospIN%frac_out,1)  /= cospIN%Npoints .OR. &
+           size(cospIN%tau_067,1)   /= cospIN%Npoints .OR. &
+           size(cospIN%emiss_11,1)  /= cospIN%Npoints .OR. &
+           size(cospgridIN%skt)     /= cospIN%Npoints .OR. &
+           size(cospgridIN%qv,1)    /= cospIN%Npoints .OR. &
+           size(cospgridIN%at,1)    /= cospIN%Npoints .OR. &
+           size(cospgridIN%phalf,1) /= cospIN%Npoints .OR. &
+           size(cospgridIN%sunlit)  /= cospIN%Npoints .OR. &
+           size(cospgridIN%pfull,1) /= cospIN%Npoints) then
           Lisccp_subcolumn = .false.
           Lisccp_column    = .false.
@@ -3403,7 +3403,7 @@
           errorMessage(nError) = 'ERROR(isccp_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%frac_out,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tau_067,2)  .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%emiss_11,2) .ne. cospIN%Ncolumns) then
+       if (size(cospIN%frac_out,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%tau_067,2)  /= cospIN%Ncolumns .OR. &
+           size(cospIN%emiss_11,2) /= cospIN%Ncolumns) then
           Lisccp_subcolumn = .false.
           Lisccp_column    = .false.
@@ -3411,11 +3411,11 @@
           errorMessage(nError) = 'ERROR(isccp_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%frac_out,3)  .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tau_067,3)   .ne. cospIN%Nlevels .OR. &
-           size(cospIN%emiss_11,3)  .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%qv,2)    .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%at,2)    .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%pfull,2) .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%phalf,2) .ne. cospIN%Nlevels+1) then
+       if (size(cospIN%frac_out,3)  /= cospIN%Nlevels .OR. &
+           size(cospIN%tau_067,3)   /= cospIN%Nlevels .OR. &
+           size(cospIN%emiss_11,3)  /= cospIN%Nlevels .OR. &
+           size(cospgridIN%qv,2)    /= cospIN%Nlevels .OR. &
+           size(cospgridIN%at,2)    /= cospIN%Nlevels .OR. &
+           size(cospgridIN%pfull,2) /= cospIN%Nlevels .OR. &
+           size(cospgridIN%phalf,2) /= cospIN%Nlevels+1) then
           Lisccp_subcolumn = .false.
           Lisccp_column    = .false.
@@ -3427,8 +3427,8 @@
     ! MISR
     if (Lmisr_subcolumn .or. Lmisr_column) then
-       if (size(cospIN%tau_067,1)        .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%sunlit)       .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%hgt_matrix,1) .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%at,1)         .ne. cospIN%Npoints) then
+       if (size(cospIN%tau_067,1)        /= cospIN%Npoints .OR. &
+           size(cospgridIN%sunlit)       /= cospIN%Npoints .OR. &
+           size(cospgridIN%hgt_matrix,1) /= cospIN%Npoints .OR. &
+           size(cospgridIN%at,1)         /= cospIN%Npoints) then
           Lmisr_subcolumn = .false.
           Lmisr_column    = .false.
@@ -3436,5 +3436,5 @@
           errorMessage(nError) = 'ERROR(misr_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%tau_067,2) .ne. cospIN%Ncolumns) then
+       if (size(cospIN%tau_067,2) /= cospIN%Ncolumns) then
           Lmisr_subcolumn = .false.
           Lmisr_column    = .false.
@@ -3442,7 +3442,7 @@
           errorMessage(nError) = 'ERROR(misr_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%tau_067,3)        .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%hgt_matrix,2) .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%at,2)         .ne. cospIN%Nlevels) then
+       if (size(cospIN%tau_067,3)        /= cospIN%Nlevels .OR. &
+           size(cospgridIN%hgt_matrix,2) /= cospIN%Nlevels .OR. &
+           size(cospgridIN%at,2)         /= cospIN%Nlevels) then
           Lmisr_subcolumn = .false.
           Lmisr_column    = .false.
@@ -3454,8 +3454,8 @@
     ! MODIS
     if (Lmodis_subcolumn .or. Lmodis_column) then
-       if (size(cospIN%fracLiq,1) .ne. cospIN%Npoints .OR. &
-           size(cospIN%tau_067,1) .ne. cospIN%Npoints .OR. &
-           size(cospIN%asym,1)    .ne. cospIN%Npoints .OR. &
-           size(cospIN%ss_alb,1)  .ne. cospIN%Npoints) then
+       if (size(cospIN%fracLiq,1) /= cospIN%Npoints .OR. &
+           size(cospIN%tau_067,1) /= cospIN%Npoints .OR. &
+           size(cospIN%asym,1)    /= cospIN%Npoints .OR. &
+           size(cospIN%ss_alb,1)  /= cospIN%Npoints) then
           Lmodis_subcolumn = .false.
           Lmodis_column    = .false.
@@ -3463,8 +3463,8 @@
           errorMessage(nError) = 'ERROR(modis_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%fracLiq,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tau_067,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%asym,2)    .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%ss_alb,2)  .ne. cospIN%Ncolumns) then
+       if (size(cospIN%fracLiq,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%tau_067,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%asym,2)    /= cospIN%Ncolumns .OR. &
+           size(cospIN%ss_alb,2)  /= cospIN%Ncolumns) then
           Lmodis_subcolumn = .false.
           Lmodis_column    = .false.
@@ -3472,8 +3472,8 @@
           errorMessage(nError) = 'ERROR(modis_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%fracLiq,3) .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tau_067,3) .ne. cospIN%Nlevels .OR. &
-           size(cospIN%asym,3)    .ne. cospIN%Nlevels .OR. &
-           size(cospIN%ss_alb,3)  .ne. cospIN%Nlevels) then
+       if (size(cospIN%fracLiq,3) /= cospIN%Nlevels .OR. &
+           size(cospIN%tau_067,3) /= cospIN%Nlevels .OR. &
+           size(cospIN%asym,3)    /= cospIN%Nlevels .OR. &
+           size(cospIN%ss_alb,3)  /= cospIN%Nlevels) then
           Lmodis_subcolumn = .false.
           Lmodis_column    = .false.
@@ -3485,8 +3485,8 @@
     ! CLOUDSAT
     if (Lcloudsat_subcolumn .or. Lcloudsat_column) then
-       if (size(cospIN%z_vol_cloudsat,1)   .ne. cospIN%Npoints .OR. &
-           size(cospIN%kr_vol_cloudsat,1)  .ne. cospIN%Npoints .OR. &
-           size(cospIN%g_vol_cloudsat,1)   .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%hgt_matrix,1)   .ne. cospIN%Npoints) then
+       if (size(cospIN%z_vol_cloudsat,1)   /= cospIN%Npoints .OR. &
+           size(cospIN%kr_vol_cloudsat,1)  /= cospIN%Npoints .OR. &
+           size(cospIN%g_vol_cloudsat,1)   /= cospIN%Npoints .OR. &
+           size(cospgridIN%hgt_matrix,1)   /= cospIN%Npoints) then
           Lcloudsat_subcolumn = .false.
           Lcloudsat_column    = .false.
@@ -3494,7 +3494,7 @@
           errorMessage(nError) = 'ERROR(cloudsat_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%z_vol_cloudsat,2)  .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%kr_vol_cloudsat,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%g_vol_cloudsat,2)  .ne. cospIN%Ncolumns) then
+       if (size(cospIN%z_vol_cloudsat,2)  /= cospIN%Ncolumns .OR. &
+           size(cospIN%kr_vol_cloudsat,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%g_vol_cloudsat,2)  /= cospIN%Ncolumns) then
           Lcloudsat_subcolumn = .false.
           Lcloudsat_column    = .false.
@@ -3502,8 +3502,8 @@
           errorMessage(nError) = 'ERROR(cloudsat_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%z_vol_cloudsat,3)  .ne. cospIN%Nlevels .OR. &
-           size(cospIN%kr_vol_cloudsat,3) .ne. cospIN%Nlevels .OR. &
-           size(cospIN%g_vol_cloudsat,3)  .ne. cospIN%Nlevels .OR. &
-           size(cospgridIN%hgt_matrix,2)  .ne. cospIN%Nlevels) then
+       if (size(cospIN%z_vol_cloudsat,3)  /= cospIN%Nlevels .OR. &
+           size(cospIN%kr_vol_cloudsat,3) /= cospIN%Nlevels .OR. &
+           size(cospIN%g_vol_cloudsat,3)  /= cospIN%Nlevels .OR. &
+           size(cospgridIN%hgt_matrix,2)  /= cospIN%Nlevels) then
           Lcloudsat_subcolumn = .false.
           Lcloudsat_column    = .false.
@@ -3515,8 +3515,8 @@
     ! GROUND LIDAR @ 532nm
     if (LgrLidar532_subcolumn .or. LgrLidar532_column) then
-       if (size(cospIN%beta_mol_grLidar532,1)    .ne. cospIN%Npoints .OR. & 
-           size(cospIN%betatot_grLidar532,1)     .ne. cospIN%Npoints .OR. &
-           size(cospIN%tau_mol_grLidar532,1)     .ne. cospIN%Npoints .OR. &
-           size(cospIN%tautot_grLidar532,1)      .ne. cospIN%Npoints) then
+       if (size(cospIN%beta_mol_grLidar532,1)    /= cospIN%Npoints .OR. &
+           size(cospIN%betatot_grLidar532,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%tau_mol_grLidar532,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_grLidar532,1)      /= cospIN%Npoints) then
           LgrLidar532_subcolumn = .false.
           LgrLidar532_column    = .false.
@@ -3524,6 +3524,6 @@
           errorMessage(nError) = 'ERROR(grLidar532_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%betatot_grLidar532,2)    .ne. cospIN%Ncolumns .OR. & 
-           size(cospIN%tautot_grLidar532,2)     .ne. cospIN%Ncolumns) then 
+       if (size(cospIN%betatot_grLidar532,2)    /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_grLidar532,2)     /= cospIN%Ncolumns) then
           LgrLidar532_subcolumn = .false. 
           LgrLidar532_column    = .false. 
@@ -3531,8 +3531,8 @@
           errorMessage(nError) = 'ERROR(grLidar532_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%beta_mol_grLidar532,2)    .ne. cospIN%Nlevels .OR. &
-           size(cospIN%betatot_grLidar532,3)     .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tau_mol_grLidar532,2)     .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tautot_grLidar532,3)      .ne. cospIN%Nlevels) then
+       if (size(cospIN%beta_mol_grLidar532,2)    /= cospIN%Nlevels .OR. &
+           size(cospIN%betatot_grLidar532,3)     /= cospIN%Nlevels .OR. &
+           size(cospIN%tau_mol_grLidar532,2)     /= cospIN%Nlevels .OR. &
+           size(cospIN%tautot_grLidar532,3)      /= cospIN%Nlevels) then
           LgrLidar532_subcolumn = .false. 
           LgrLidar532_column    = .false. 
@@ -3544,8 +3544,8 @@
     ! ATLID
     if (Latlid_subcolumn .or. Latlid_column) then
-       if (size(cospIN%beta_mol_atlid,1)    .ne. cospIN%Npoints .OR. &
-           size(cospIN%betatot_atlid,1)     .ne. cospIN%Npoints .OR. &
-           size(cospIN%tau_mol_atlid,1)     .ne. cospIN%Npoints .OR. & 
-           size(cospIN%tautot_atlid,1)      .ne. cospIN%Npoints) then 
+       if (size(cospIN%beta_mol_atlid,1)    /= cospIN%Npoints .OR. &
+           size(cospIN%betatot_atlid,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%tau_mol_atlid,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_atlid,1)      /= cospIN%Npoints) then
           Latlid_subcolumn = .false. 
           Latlid_column    = .false. 
@@ -3553,6 +3553,6 @@
           errorMessage(nError) = 'ERROR(atlid_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%betatot_atlid,2)    .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tautot_atlid,2)     .ne. cospIN%Ncolumns) then 
+       if (size(cospIN%betatot_atlid,2)    /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_atlid,2)     /= cospIN%Ncolumns) then
           Latlid_subcolumn = .false.
           Latlid_column    = .false.
@@ -3560,8 +3560,8 @@
           errorMessage(nError) = 'ERROR(atlid_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%beta_mol_atlid,2)    .ne. cospIN%Nlevels .OR. &
-           size(cospIN%betatot_atlid,3)     .ne. cospIN%Nlevels .OR. & 
-           size(cospIN%tau_mol_atlid,2)     .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tautot_atlid,3)      .ne. cospIN%Nlevels) then 
+       if (size(cospIN%beta_mol_atlid,2)    /= cospIN%Nlevels .OR. &
+           size(cospIN%betatot_atlid,3)     /= cospIN%Nlevels .OR. &
+           size(cospIN%tau_mol_atlid,2)     /= cospIN%Nlevels .OR. &
+           size(cospIN%tautot_atlid,3)      /= cospIN%Nlevels) then
           Latlid_subcolumn = .false.
           Latlid_column    = .false. 
@@ -3573,12 +3573,12 @@
     ! CALIPSO
     if (Lcalipso_subcolumn .or. Lcalipso_column) then
-       if (size(cospIN%beta_mol_calipso,1)    .ne. cospIN%Npoints .OR. &
-           size(cospIN%betatot_calipso,1)     .ne. cospIN%Npoints .OR. &
-           size(cospIN%betatot_liq_calipso,1) .ne. cospIN%Npoints .OR. &
-           size(cospIN%betatot_ice_calipso,1) .ne. cospIN%Npoints .OR. &
-           size(cospIN%tau_mol_calipso,1)     .ne. cospIN%Npoints .OR. &
-           size(cospIN%tautot_calipso,1)      .ne. cospIN%Npoints .OR. &
-           size(cospIN%tautot_liq_calipso,1)  .ne. cospIN%Npoints .OR. &
-           size(cospIN%tautot_ice_calipso,1)  .ne. cospIN%Npoints) then
+       if (size(cospIN%beta_mol_calipso,1)    /= cospIN%Npoints .OR. &
+           size(cospIN%betatot_calipso,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%betatot_liq_calipso,1) /= cospIN%Npoints .OR. &
+           size(cospIN%betatot_ice_calipso,1) /= cospIN%Npoints .OR. &
+           size(cospIN%tau_mol_calipso,1)     /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_calipso,1)      /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_liq_calipso,1)  /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_ice_calipso,1)  /= cospIN%Npoints) then
           Lcalipso_subcolumn = .false.
           Lcalipso_column    = .false.
@@ -3586,10 +3586,10 @@
           errorMessage(nError) = 'ERROR(calipso_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%betatot_calipso,2)     .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%betatot_liq_calipso,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%betatot_ice_calipso,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tautot_calipso,2)      .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tautot_liq_calipso,2)  .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tautot_ice_calipso,2)  .ne. cospIN%Ncolumns) then
+       if (size(cospIN%betatot_calipso,2)     /= cospIN%Ncolumns .OR. &
+           size(cospIN%betatot_liq_calipso,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%betatot_ice_calipso,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_calipso,2)      /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_liq_calipso,2)  /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_ice_calipso,2)  /= cospIN%Ncolumns) then
           Lcalipso_subcolumn = .false.
           Lcalipso_column    = .false.
@@ -3597,12 +3597,12 @@
           errorMessage(nError) = 'ERROR(calipso_simulator): The number of sub-columns in the input fields are inconsistent'
        endif
-       if (size(cospIN%beta_mol_calipso,2)    .ne. cospIN%Nlevels .OR. &
-           size(cospIN%betatot_calipso,3)     .ne. cospIN%Nlevels .OR. &
-           size(cospIN%betatot_liq_calipso,3) .ne. cospIN%Nlevels .OR. &
-           size(cospIN%betatot_ice_calipso,3) .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tau_mol_calipso,2)     .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tautot_calipso,3)      .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tautot_liq_calipso,3)  .ne. cospIN%Nlevels .OR. &
-           size(cospIN%tautot_ice_calipso,3)  .ne. cospIN%Nlevels) then
+       if (size(cospIN%beta_mol_calipso,2)    /= cospIN%Nlevels .OR. &
+           size(cospIN%betatot_calipso,3)     /= cospIN%Nlevels .OR. &
+           size(cospIN%betatot_liq_calipso,3) /= cospIN%Nlevels .OR. &
+           size(cospIN%betatot_ice_calipso,3) /= cospIN%Nlevels .OR. &
+           size(cospIN%tau_mol_calipso,2)     /= cospIN%Nlevels .OR. &
+           size(cospIN%tautot_calipso,3)      /= cospIN%Nlevels .OR. &
+           size(cospIN%tautot_liq_calipso,3)  /= cospIN%Nlevels .OR. &
+           size(cospIN%tautot_ice_calipso,3)  /= cospIN%Nlevels) then
           Lcalipso_subcolumn = .false.
           Lcalipso_column    = .false.
@@ -3614,6 +3614,6 @@
     ! PARASOL
     if (Lparasol_subcolumn .or. Lparasol_column) then
-       if (size(cospIN%tautot_S_liq,1) .ne. cospIN%Npoints .OR. &
-           size(cospIN%tautot_S_ice,1) .ne. cospIN%Npoints) then
+       if (size(cospIN%tautot_S_liq,1) /= cospIN%Npoints .OR. &
+           size(cospIN%tautot_S_ice,1) /= cospIN%Npoints) then
           Lparasol_subcolumn = .false.
           Lparasol_column    = .false.
@@ -3621,6 +3621,6 @@
           errorMessage(nError) = 'ERROR(parasol_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospIN%tautot_S_liq,2) .ne. cospIN%Ncolumns .OR. &
-           size(cospIN%tautot_S_ice,2) .ne. cospIN%Ncolumns) then
+       if (size(cospIN%tautot_S_liq,2) /= cospIN%Ncolumns .OR. &
+           size(cospIN%tautot_S_ice,2) /= cospIN%Ncolumns) then
           Lparasol_subcolumn = .false.
           Lparasol_column    = .false.
@@ -3632,25 +3632,25 @@
     ! RTTOV
     if (Lrttov_column) then
-       if (size(cospgridIN%pfull,1)           .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%at,1)              .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%qv,1)              .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%hgt_matrix_half,1) .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%u_sfc)             .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%v_sfc)             .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%skt)               .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%phalf,1)           .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%qv,1)              .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%land)              .ne. cospIN%Npoints .OR. &
-           size(cospgridIN%lat)               .ne. cospIN%Npoints) then
+       if (size(cospgridIN%pfull,1)           /= cospIN%Npoints .OR. &
+           size(cospgridIN%at,1)              /= cospIN%Npoints .OR. &
+           size(cospgridIN%qv,1)              /= cospIN%Npoints .OR. &
+           size(cospgridIN%hgt_matrix_half,1) /= cospIN%Npoints .OR. &
+           size(cospgridIN%u_sfc)             /= cospIN%Npoints .OR. &
+           size(cospgridIN%v_sfc)             /= cospIN%Npoints .OR. &
+           size(cospgridIN%skt)               /= cospIN%Npoints .OR. &
+           size(cospgridIN%phalf,1)           /= cospIN%Npoints .OR. &
+           size(cospgridIN%qv,1)              /= cospIN%Npoints .OR. &
+           size(cospgridIN%land)              /= cospIN%Npoints .OR. &
+           size(cospgridIN%lat)               /= cospIN%Npoints) then
           Lrttov_column    = .false.
           nError=nError+1
           errorMessage(nError) = 'ERROR(rttov_simulator): The number of points in the input fields are inconsistent'
        endif
-       if (size(cospgridIN%pfull,2)           .ne. cospIN%Nlevels   .OR. &
-           size(cospgridIN%at,2)              .ne. cospIN%Nlevels   .OR. &
-           size(cospgridIN%qv,2)              .ne. cospIN%Nlevels   .OR. &
-           size(cospgridIN%hgt_matrix_half,2) .ne. cospIN%Nlevels+1 .OR. &
-           size(cospgridIN%phalf,2)           .ne. cospIN%Nlevels+1 .OR. &
-           size(cospgridIN%qv,2)              .ne. cospIN%Nlevels) then
+       if (size(cospgridIN%pfull,2)           /= cospIN%Nlevels   .OR. &
+           size(cospgridIN%at,2)              /= cospIN%Nlevels   .OR. &
+           size(cospgridIN%qv,2)              /= cospIN%Nlevels   .OR. &
+           size(cospgridIN%hgt_matrix_half,2) /= cospIN%Nlevels+1 .OR. &
+           size(cospgridIN%phalf,2)           /= cospIN%Nlevels+1 .OR. &
+           size(cospgridIN%qv,2)              /= cospIN%Nlevels) then
           Lrttov_column    = .false.
           nError=nError+1
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_stats.F90	(revision 5082)
@@ -210,5 +210,5 @@
           ! 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) == R_GROUND)
              j = j+1
           enddo
@@ -217,7 +217,7 @@
           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
-             if (Ze_tot(pr,i,j) .lt. -30.) then  !radar can't detect cloud
-                if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then  !lidar sense cloud
+             if ((sc_ratio <= s_att) .and. (flag_sat == 0)) flag_sat = j
+             if (Ze_tot(pr,i,j) < -30.) then  !radar can't detect cloud
+                if ( (sc_ratio > s_cld) .or. (flag_sat == j) ) then  !lidar sense cloud
                    lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
                    flag_cld=1
@@ -226,10 +226,10 @@
                 flag_cld=1
                 flag_radarcld=1
-                if (j .gt. j_1km) flag_radarcld_no1km=1              
+                if (j > j_1km) flag_radarcld_no1km=1
              endif
           enddo !levels
-          if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1._wp
-          if (flag_radarcld .eq. 1) radar_tcc(pr)=radar_tcc(pr)+1.
-          if (flag_radarcld_no1km .eq. 1) radar_tcc2(pr)=radar_tcc2(pr)+1.        
+          if (flag_cld == 1) tcc(pr)=tcc(pr)+1._wp
+          if (flag_radarcld == 1) radar_tcc(pr)=radar_tcc(pr)+1.
+          if (flag_radarcld_no1km == 1) radar_tcc2(pr)=radar_tcc2(pr)+1.
        enddo !columns
     enddo !points
@@ -267,6 +267,6 @@
     
     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
+       hist1D(ij-1) = count(var >= bins(ij-1) .and. var < bins(ij))
+       if (count(var == R_GROUND) >= 1) hist1D(ij-1)=R_UNDEF
     enddo
     
@@ -300,6 +300,6 @@
     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))        
+          jointHist(ij-1,ik-1)=count(var1 >= bin1(ij-1) .and. var1 < bin1(ij) .and. &
+               var2 >= bin2(ik-1) .and. var2 < bin2(ik))
        enddo
     enddo
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/icarus.F90	(revision 5082)
@@ -134,5 +134,5 @@
     ! ##########################################################################
     
-    if (debugcol.ne.0) then
+    if (debugcol/=0) then
        do j=1,npoints,debugcol
           
@@ -140,5 +140,5 @@
           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
+             where(levmatch(j,1:ncol) == ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
           enddo
           
@@ -224,5 +224,5 @@
 
     ! Set tropopause values
-    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then 
+    if (isccp_top_height == 1 .or. isccp_top_height == 3) then
        ptrop(1:npoints)     = 5000._wp
        attropmin(1:npoints) = 400._wp
@@ -232,7 +232,7 @@
 
        do ilev=1,nlev
-          where(pfull(1:npoints,ilev) .lt. 40000. .and. &
-                pfull(1:npoints,ilev) .gt.  5000. .and. &
-                at(1:npoints,ilev)    .lt. attropmin(1:npoints))
+          where(pfull(1:npoints,ilev) < 40000. .and. &
+                pfull(1:npoints,ilev) >  5000. .and. &
+                at(1:npoints,ilev)    < attropmin(1:npoints))
              ptrop(1:npoints)     = pfull(1:npoints,ilev)
              attropmin(1:npoints) = at(1:npoints,ilev)
@@ -244,9 +244,9 @@
        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))
+               at(1:npoints,ilev) > atmax(1:npoints) .and. ilev  >= itrop(1:npoints))
        enddo
     end if
   
-    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
+    if (isccp_top_height == 1 .or. isccp_top_height == 3) then
        ! ############################################################################
        !                        Clear-sky radiance calculation
@@ -308,5 +308,5 @@
              dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
                                          1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
-                                         demIN(1:npoints,ibox,ilev) .eq. 0)
+                                         demIN(1:npoints,ibox,ilev) == 0)
 
              ! Increase TOA flux emitted from layer
@@ -348,7 +348,7 @@
           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
+          if (isccp_top_height == 1) then
              do j=1,npoints  
-                if (transmax(j) .gt. 0.001 .and.  transmax(j) .le. 0.9999999) then
+                if (transmax(j) > 0.001 .and.  transmax(j) <= 0.9999999) then
                    fluxtopinit(j) = fluxtop(j,ibox)
                    tauir(j) = tau(j,ibox)/2.13_wp
@@ -357,11 +357,11 @@
              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
+                   if (tau(j,ibox) > (tauchk)) then
+                      if (transmax(j) > 0.001 .and.  transmax(j) <= 0.9999999) then
                          emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j)  )
                          fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
                          fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
                          tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
-                         if (tb(j,ibox) .gt. 260.) then
+                         if (tb(j,ibox) > 260.) then
                             tauir(j) = tau(j,ibox) / 2.56_wp
                          end if
@@ -373,7 +373,7 @@
 
           ! Cloud-top temperature
-          where(tau(1:npoints,ibox) .gt. tauchk)
+          where(tau(1:npoints,ibox) > tauchk)
              tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox))))
-             where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints))
+             where (isccp_top_height == 1 .and. tauir(1:npoints) < taumin(1:npoints))
                 tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp 
                 tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
@@ -382,5 +382,5 @@
           
           ! Clear-sky brightness temperature
-          where(tau(1:npoints,ibox) .le. tauchk) 
+          where(tau(1:npoints,ibox) <= tauchk)
              tb(1:npoints,ibox) = meantbclr(1:npoints)
           endwhere
@@ -399,16 +399,16 @@
     do ibox=1,ncol
        !segregate according to optical thickness
-       if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then  
+       if (isccp_top_height == 1 .or. isccp_top_height == 3) then
           
           ! Find level whose temperature most closely matches brightness temperature
           nmatch(1:npoints)=0
           do k1=1,nlev-1
-             ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)        
+             ilev = merge(nlev-k1,k1,isccp_top_height_direction == 2)
              do j=1,npoints 
-                if (ilev           .ge. itrop(j)     .and. &
-                     ((at(j,ilev)  .ge. tb(j,ibox)   .and. &  
-                      at(j,ilev+1) .le. tb(j,ibox))  .or.  &
-                      (at(j,ilev)  .le. tb(j,ibox)   .and. &
-                      at(j,ilev+1) .ge. tb(j,ibox)))) then 
+                if (ilev           >= itrop(j)     .and. &
+                     ((at(j,ilev)  >= tb(j,ibox)   .and. &
+                      at(j,ilev+1) <= tb(j,ibox))  .or.  &
+                      (at(j,ilev)  <= tb(j,ibox)   .and. &
+                      at(j,ilev+1) >= tb(j,ibox)))) then
                    nmatch(j)=nmatch(j)+1
                    match(j,nmatch(j))=ilev
@@ -418,5 +418,5 @@
 
           do j=1,npoints 
-             if (nmatch(j) .ge. 1) then
+             if (nmatch(j) >= 1) then
                 k1 = match(j,nmatch(j))
                 k2 = k1 + 1
@@ -426,11 +426,11 @@
                 logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
                 ptop(j,ibox) = exp(logp)
-                levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox)))
+                levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) < abs(pfull(j,k2)-ptop(j,ibox)))
              else
-                if (tb(j,ibox) .le. attrop(j)) then
+                if (tb(j,ibox) <= attrop(j)) then
                    ptop(j,ibox)=ptrop(j)
                    levmatch(j,ibox)=itrop(j)
                 end if
-                if (tb(j,ibox) .ge. atmax(j)) then
+                if (tb(j,ibox) >= atmax(j)) then
                    ptop(j,ibox)=pfull(j,nlev)
                    levmatch(j,ibox)=nlev
@@ -441,5 +441,5 @@
           ptop(1:npoints,ibox)=0.
           do ilev=1,nlev
-             where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0))
+             where((ptop(1:npoints,ibox) == 0. ) .and.(frac_out(1:npoints,ibox,ilev) /= 0))
                 ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
                 levmatch(1:npoints,ibox)=ilev
@@ -447,5 +447,5 @@
           end do
        end if
-       where(tau(1:npoints,ibox) .le. tauchk)
+       where(tau(1:npoints,ibox) <= tauchk)
           ptop(1:npoints,ibox)=0._wp
           levmatch(1:npoints,ibox)=0._wp
@@ -460,6 +460,6 @@
     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
+          if (tau(j,ibox) > (tauchk) .and. ptop(j,ibox) > 0.) then
+             if (sunlit(j)==1 .or. isccp_top_height == 3) then
                 boxtau(j,ibox) = tau(j,ibox)
                 boxptop(j,ibox) = ptop(j,ibox)!/100._wp
@@ -508,5 +508,5 @@
     !                           Brightness Temperature
     ! ####################################################################################
-    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
+    if (isccp_top_height == 1 .or. isccp_top_height == 3) then
        meantb(1:npoints)=sum(boxttop,2)/ncol
     else
@@ -535,5 +535,5 @@
        do ilev2=1,7
           do j=1,npoints ! 
-             if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then 
+             if (sunlit(j)==1 .or. isccp_top_height == 3) then
                 fq_isccp(j,ilev,ilev2)= 0.
 	     else 
@@ -546,5 +546,5 @@
     
     ! Reset variables need for averaging cloud properties
-    where(sunlit .eq. 1 .or. isccp_top_height .eq. 3)
+    where(sunlit == 1 .or. isccp_top_height == 3)
        totalcldarea(1:npoints)  = 0._wp
        meanalbedocld(1:npoints) = 0._wp
@@ -561,8 +561,8 @@
     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.)
+       box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) > tauchk .and. boxptop(j,1:ncol) > 0.)
 
        ! Compute joint histogram and column quantities for points that are sunlit and cloudy
-       if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then 
+       if (sunlit(j) ==1 .or. isccp_top_height == 3) then
           ! Joint-histogram
           call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, &
@@ -572,10 +572,10 @@
           
           ! Column cloud area
-          totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol
+          totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin))/ncol
              
           ! Subcolumn cloud albedo
           !albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),&
           !     0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
-          where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
+          where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin)
              albedocld(j,1:ncol) = (boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp)
           elsewhere
@@ -587,10 +587,10 @@
           
           ! Column cloud top pressure
-          meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol
+          meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) > isccp_taumin)/ncol
        endif
     enddo
     
     ! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
-    where(totalcldarea(1:npoints) .gt. 0)
+    where(totalcldarea(1:npoints) > 0)
        meanptop(1:npoints)      = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
        meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
@@ -609,6 +609,6 @@
 
     ! Represent in percent
-    where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp
-    where(fq_isccp     .ne. output_missing_value) fq_isccp     = fq_isccp*100._wp
+    where(totalcldarea /= output_missing_value) totalcldarea = totalcldarea*100._wp
+    where(fq_isccp     /= output_missing_value) fq_isccp     = fq_isccp*100._wp
     
     
@@ -634,8 +634,8 @@
     varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
     do j=1,dim2
-       where(flag(:,j,:) .eq. 1)
+       where(flag(:,j,:) == 1)
           varOUT(:,j,:) = varIN2
        endwhere
-       where(flag(:,j,:) .eq. 2)
+       where(flag(:,j,:) == 2)
           varOUT(:,j,:) = varIN1
        endwhere
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lidar_simulator.F90	(revision 5082)
@@ -231,5 +231,5 @@
           ! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
           ! Upper layer
-          WHERE(tautot(1:npoints,icol,1) .gt. 0)
+          WHERE(tautot(1:npoints,icol,1) > 0)
              pnorm_perp_tot(1:npoints,icol,1) = (beta_perp_ice(1:npoints,icol,1)+           &
                   beta_perp_liq(1:npoints,icol,1)-                                          &
@@ -251,8 +251,8 @@
              ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following 
              ! equations:
-             WHERE (pnorm(1:npoints,icol,k) .eq. 0)
+             WHERE (pnorm(1:npoints,icol,k) == 0)
                 pnorm_perp_tot(1:npoints,icol,k)=0._wp
              ELSEWHERE
-                where(tautot_lay(1:npoints) .gt. 0.)
+                where(tautot_lay(1:npoints) > 0.)
                    pnorm_perp_tot(1:npoints,icol,k) = (beta_perp_ice(1:npoints,icol,k)+     &
                         beta_perp_liq(1:npoints,icol,k)-(beta_mol(1:npoints,k)/(1._wp+1._wp/  &
@@ -358,7 +358,7 @@
     latlid = .false.
     lgrlidar532 = .false.
-    if (platform .eq. 'calipso') lcalipso=.true.
-    if (platform .eq. 'atlid') latlid=.true.
-    if (platform .eq. 'grlidar532') lgrlidar532=.true.
+    if (platform == 'calipso') lcalipso=.true.
+    if (platform == 'atlid') latlid=.true.
+    if (platform == 'grlidar532') lgrlidar532=.true.
         
     ! Vertically regrid input data
@@ -400,6 +400,6 @@
        do ic = 1, ncol
           pnorm_c = pnormFlip(:,ic,:)
-          where ((pnorm_c .lt. xmax) .and. (betamolFlip(:,1,:) .lt. xmax) .and.          &
-                (betamolFlip(:,1,:) .gt. 0.0 ))
+          where ((pnorm_c < xmax) .and. (betamolFlip(:,1,:) < xmax) .and.          &
+                (betamolFlip(:,1,:) > 0.0 ))
              x3d_c = pnorm_c/betamolFlip(:,1,:)
           elsewhere
@@ -429,5 +429,5 @@
        do ic = 1, ncol
           pnorm_c = pnorm(:,ic,:)
-          where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
+          where ((pnorm_c<xmax) .and. (pmol<xmax) .and. (pmol> 0.0 ))
              x3d_c = pnorm_c/pmol
           elsewhere
@@ -463,5 +463,5 @@
           enddo
        enddo
-       where(cfad2 .ne. R_UNDEF) cfad2=cfad2/ncol
+       where(cfad2 /= R_UNDEF) cfad2=cfad2/ncol
     endif 
     
@@ -501,5 +501,5 @@
     do k=2,nlev
        tautot_lay(:) = tau(:,k)-tau(:,k-1) 
-       WHERE (tautot_lay(:) .gt. 0.)
+       WHERE (tautot_lay(:) > 0.)
           pnorm(:,k) = beta(:,k)*EXP(-2._wp*tau(:,k-1)) /&
                (2._wp*tautot_lay(:))*(1._wp-EXP(-2._wp*tautot_lay(:)))
@@ -529,6 +529,6 @@
     do k=2,nlev
        tautot_lay(:) = tau(:,k)-tau(:,k-1)       
-       WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. epsrealwp )
-          WHERE (tautot_lay(:) .gt. 0.)
+       WHERE ( EXP(-2._wp*tau(:,k-1)) > epsrealwp )
+          WHERE (tautot_lay(:) > 0.)
              beta(:,k) = pnorm(:,k)/ EXP(-2._wp*tau(:,k-1))* &
                   (2._wp*tautot_lay(:))/(1._wp-exp(-2._wp*tautot_lay(:)))
@@ -650,5 +650,5 @@
     do k=1,Nlevels
        ! Cloud detection at subgrid-scale:
-       where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
+       where ((x(:,:,k) > S_cld) .and. (x(:,:,k) /= undef) )
           cldy(:,:,k)=1._wp
        elsewhere
@@ -657,5 +657,5 @@
        
        ! Number of usefull sub-columns:
-       where ((x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
+       where ((x(:,:,k) > S_att) .and. (x(:,:,k) /= undef) )
           srok(:,:,k)=1._wp
        elsewhere
@@ -677,7 +677,7 @@
              ! 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
+             if(srok(ip,ic,k)>0.)then
                 do itemp=1,Ntemp
-                   if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
+                   if( (tmp(ip,k)>=tempmod(itemp)).and.(tmp(ip,k)<tempmod(itemp+1)) )then
                       lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1._wp
                    endif
@@ -685,7 +685,7 @@
              endif
              
-             if(cldy(ip,ic,k).eq.1.)then
+             if(cldy(ip,ic,k)==1.)then
                 do itemp=1,Ntemp 
-                   if( (tmp(ip,k) .ge. tempmod(itemp)).and.(tmp(ip,k) .lt. tempmod(itemp+1)) )then
+                   if( (tmp(ip,k) >= tempmod(itemp)).and.(tmp(ip,k) < tempmod(itemp+1)) )then
                       lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1._wp
                    endif
@@ -695,7 +695,7 @@
              iz=1
              p1 = pplay(ip,k)
-             if ( p1.gt.0. .and. p1.lt.(440._wp*100._wp)) then ! high clouds
+             if ( p1>0. .and. p1<(440._wp*100._wp)) then ! high clouds
                 iz=3
-             else if(p1.ge.(440._wp*100._wp) .and. p1.lt.(680._wp*100._wp)) then ! mid clouds
+             else if(p1>=(440._wp*100._wp) .and. p1<(680._wp*100._wp)) then ! mid clouds
                 iz=2
              endif
@@ -714,5 +714,5 @@
     
     ! Grid-box 3D cloud fraction
-    where ( nsub(:,:).gt.0.0 )
+    where ( nsub(:,:)>0.0 )
        lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
     elsewhere
@@ -729,5 +729,5 @@
        enddo
     enddo
-    where (nsublayer(:,:) .gt. 0.0)
+    where (nsublayer(:,:) > 0.0)
        cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
     elsewhere
@@ -748,5 +748,5 @@
 
              ! Avoid zero values
-             if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
+             if( (cldy(i,ncol,nlev)==1.) .and. (ATBperp(i,ncol,nlev)>0.) )then
                 ! Computation of the ATBperp along the phase discrimination line
                 ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
@@ -756,8 +756,8 @@
                 ! 4.1.a) Ice: ATBperp above the phase discrimination line
                 ! ########################################################################
-                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
+                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) >= 0.)then ! Ice clouds
 
                    ! ICE with temperature above 273,15°K = Liquid (false ice)
-                   if(tmp(i,nlev) .gt. 273.15) then ! Temperature above 273,15 K
+                   if(tmp(i,nlev) > 273.15) then ! Temperature above 273,15 K
                      ! Liquid: False ice corrected by the temperature to Liquid
                       lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! False ice detection ==> added to Liquid
@@ -767,7 +767,7 @@
                                                                               ! to classify the phase cloud
                       cldlayphase(i,ncol,4,2) = 1. ! tot cloud
-                      if (p1 .gt. 0. .and. p1.lt.(440._wp*100._wp)) then ! high cloud
+                      if (p1 > 0. .and. p1<(440._wp*100._wp)) then ! high cloud
                          cldlayphase(i,ncol,3,2) = 1._wp
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then ! mid cloud
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then ! mid cloud
                          cldlayphase(i,ncol,2,2) = 1._wp
                       else ! low cloud
@@ -776,8 +776,8 @@
                       cldlayphase(i,ncol,4,5) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then 
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,5) = 1._wp
                       ! Middle cloud
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,5) = 1._wp
                       ! Low cloud
@@ -791,8 +791,8 @@
                       cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud 
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then 
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,1) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then 
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,1) = 1._wp
                       ! Low cloud
@@ -806,13 +806,13 @@
                 else
                    ! Liquid with temperature above 231,15°K
-                   if(tmp(i,nlev) .gt. 231.15_wp) then
+                   if(tmp(i,nlev) > 231.15_wp) then
                       lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
                       tmpl(i,ncol,nlev)       = tmp(i,nlev)
                       cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,2) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,2) = 1._wp
                       ! Low cloud   
@@ -827,8 +827,8 @@
                       cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,4) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,4) = 1._wp
                       ! Low cloud
@@ -838,8 +838,8 @@
                       cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,1) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,1) = 1._wp
                       ! Low cloud   
@@ -859,5 +859,5 @@
              p1 = pplay(i,nlev)
 
-             if((cldy(i,ncol,nlev) .eq. 1.) .and. (ATBperp(i,ncol,nlev) .gt. 0.) )then
+             if((cldy(i,ncol,nlev) == 1.) .and. (ATBperp(i,ncol,nlev) > 0.) )then
                 ! Computation of the ATBperp of the phase discrimination line
                 ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
@@ -868,6 +868,6 @@
                 ! ########################################################################
                 ! ICE with temperature above 273,15°K = Liquid (false ice)
-                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
-                   if(tmp(i,nlev) .gt. 273.15)then
+                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) >= 0.)then ! Ice clouds
+                   if(tmp(i,nlev) > 273.15)then
                       lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! false ice ==> liq
                       tmpl(i,ncol,nlev)       = tmp(i,nlev)
@@ -875,8 +875,8 @@
                       cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then 
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,2) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then 
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,2) = 1._wp
                       ! Low cloud
@@ -887,8 +887,8 @@
                       cldlayphase(i,ncol,4,5) = 1. ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,5) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,5) = 1._wp
                       ! Low cloud   
@@ -902,8 +902,8 @@
                       cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,1) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt.(680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 <(680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,1) = 1._wp
                       ! Low cloud   
@@ -918,13 +918,13 @@
                 else
                    ! Liquid with temperature above 231,15°K
-                   if(tmp(i,nlev) .gt. 231.15)then
+                   if(tmp(i,nlev) > 231.15)then
                       lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
                       tmpl(i,ncol,nlev)       = tmp(i,nlev)
                       cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,2) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,2) = 1._wp
                       ! Low cloud   
@@ -939,8 +939,8 @@
                       cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,4) = 1._wp
                       ! Middle   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,4) = 1._wp
                       ! Low cloud   
@@ -951,8 +951,8 @@
                       cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
                       ! High cloud
-                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                      if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                          cldlayphase(i,ncol,3,1) = 1._wp
                       ! Middle cloud   
-                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                      else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                          cldlayphase(i,ncol,2,1) = 1._wp
                       ! Low cloud   
@@ -966,5 +966,5 @@
                 
                 ! Find the level of the highest cloud with SR>30
-                if(x(i,ncol,nlev) .gt. S_cld_att) then ! SR > 30.
+                if(x(i,ncol,nlev) > S_cld_att) then ! SR > 30.
                     toplvlsat = nlev+1
                     goto 99
@@ -978,16 +978,16 @@
           ! see Cesana and Chepfer 2013 Sect.III.2
           ! ##############################################################################
-          if(toplvlsat.ne.0) then
+          if(toplvlsat/=0) then
              do nlev = toplvlsat,Nlevels
                 p1 = pplay(i,nlev)
-                if(cldy(i,ncol,nlev).eq.1.)then
+                if(cldy(i,ncol,nlev)==1.)then
                    lidarcldphase(i,nlev,3) = lidarcldphase(i,nlev,3)+1._wp
                    tmpu(i,ncol,nlev)       = tmp(i,nlev)
                    cldlayphase(i,ncol,4,3) = 1._wp ! tot cloud
                    ! High cloud
-                   if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
+                   if (p1 > 0. .and. p1 < (440._wp*100._wp)) then
                       cldlayphase(i,ncol,3,3) = 1._wp
                    ! Middle cloud   
-                   else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
+                   else if(p1 >= (440._wp*100._wp) .and. p1 < (680._wp*100._wp)) then
                       cldlayphase(i,ncol,2,3) = 1._wp
                    ! Low cloud   
@@ -1008,5 +1008,5 @@
     ! Compute the Ice percentage in cloud = ice/(ice+liq) as a function of the occurrences
     lidarcldphasetmp(:,:) = lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
-    WHERE (lidarcldphasetmp(:,:) .gt. 0.)
+    WHERE (lidarcldphasetmp(:,:) > 0.)
        lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
     ELSEWHERE
@@ -1016,5 +1016,5 @@
     ! Compute Phase 3D Cloud Fraction
     !WHERE (nsub(:,Nlevels:1:-1) .gt. 0.0 )
-    WHERE (nsub(:,:) .gt. 0.0 )  
+    WHERE (nsub(:,:) > 0.0 )
        lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
        lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
@@ -1049,5 +1049,5 @@
     ! Compute the Ice percentage in cloud = ice/(ice+liq)
     cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
-    WHERE (cldlayerphasetmp(:,:).gt. 0.)
+    WHERE (cldlayerphasetmp(:,:)> 0.)
        cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
     ELSEWHERE
@@ -1056,5 +1056,5 @@
     
     do i=1,Nphase-1
-       WHERE ( cldlayerphasesum(:,:).gt.0.0 )
+       WHERE ( cldlayerphasesum(:,:)>0.0 )
           cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
        ENDWHERE
@@ -1065,10 +1065,10 @@
           checkcldlayerphase=0.
           checkcldlayerphase2=0.
-          if (cldlayerphasesum(i,iz) .gt. 0.0 )then
+          if (cldlayerphasesum(i,iz) > 0.0 )then
              do ic=1,Nphase-3
                 checkcldlayerphase = checkcldlayerphase+cldlayerphase(i,iz,ic)
              enddo
              checkcldlayerphase2 = cldlayer(i,iz)-checkcldlayerphase
-             if((checkcldlayerphase2 .gt. 0.01) .or. (checkcldlayerphase2 .lt. -0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
+             if((checkcldlayerphase2 > 0.01) .or. (checkcldlayerphase2 < -0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
           endif
        enddo
@@ -1076,5 +1076,5 @@
     
     do i=1,Nphase-1
-       WHERE (nsublayer(:,:) .eq. 0.0)
+       WHERE (nsublayer(:,:) == 0.0)
           cldlayerphase(:,:,i) = undef
        ENDWHERE
@@ -1086,14 +1086,14 @@
           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
+                if(tmpi(i,ncol,nlev)>0.)then
+                   if((tmpi(i,ncol,nlev) >= tempmod(itemp)) .and. (tmpi(i,ncol,nlev) < tempmod(itemp+1)) )then
                       lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1._wp
                    endif
-                elseif(tmpl(i,ncol,nlev) .gt. 0.)then
-                   if((tmpl(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpl(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
+                elseif(tmpl(i,ncol,nlev) > 0.)then
+                   if((tmpl(i,ncol,nlev) >= tempmod(itemp)) .and. (tmpl(i,ncol,nlev) < tempmod(itemp+1)) )then
                       lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1._wp
                    endif
-                elseif(tmpu(i,ncol,nlev) .gt. 0.)then
-                   if((tmpu(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpu(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
+                elseif(tmpu(i,ncol,nlev) > 0.)then
+                   if((tmpu(i,ncol,nlev) >= tempmod(itemp)) .and. (tmpu(i,ncol,nlev) < tempmod(itemp+1)) )then
                       lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1._wp
                    endif
@@ -1118,5 +1118,5 @@
     ! Compute the Ice percentage in cloud = ice/(ice+liq)
     sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)    
-    WHERE(sumlidarcldtemp(:,:) .gt. 0.)
+    WHERE(sumlidarcldtemp(:,:) > 0.)
        lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
     ELSEWHERE
@@ -1125,5 +1125,5 @@
     
     do i=1,4
-       WHERE(lidarcldtempind(:,:) .gt. 0.)
+       WHERE(lidarcldtempind(:,:) > 0.)
           lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
        ELSEWHERE
@@ -1193,5 +1193,5 @@
     do k=1,Nlevels
        ! Cloud detection at subgrid-scale:
-       where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
+       where ((x(:,:,k) > S_cld) .and. (x(:,:,k) /= undef) )
           cldy(:,:,k)=1._wp
        elsewhere
@@ -1200,5 +1200,5 @@
        
        ! Number of usefull sub-columns:
-       where ((x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
+       where ((x(:,:,k) > S_att) .and. (x(:,:,k) /= undef) )
           srok(:,:,k)=1._wp
        elsewhere
@@ -1216,7 +1216,7 @@
              iz=1
              p1 = pplay(ip,k)
-             if ( p1.gt.0. .and. p1.lt.(440._wp*100._wp)) then ! high clouds
+             if ( p1>0. .and. p1<(440._wp*100._wp)) then ! high clouds
                 iz=3
-             else if(p1.ge.(440._wp*100._wp) .and. p1.lt.(680._wp*100._wp)) then ! mid clouds
+             else if(p1>=(440._wp*100._wp) .and. p1<(680._wp*100._wp)) then ! mid clouds
                 iz=2
              endif
@@ -1235,5 +1235,5 @@
     
     ! Grid-box 3D cloud fraction
-    where ( nsub(:,:).gt.0.0 )
+    where ( nsub(:,:)>0.0 )
        lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
     elsewhere
@@ -1250,5 +1250,5 @@
        enddo
     enddo
-    where (nsublayer(:,:) .gt. 0.0)
+    where (nsublayer(:,:) > 0.0)
        cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
     elsewhere
@@ -1344,5 +1344,5 @@
     do k=1,Nlevels
        ! Cloud detection at subgrid-scale:
-       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
+       where ( (x(:,:,k) > S_cld) .and. (x(:,:,k) /= undef) )
           cldy(:,:,k)=1._wp
        elsewhere
@@ -1350,5 +1350,5 @@
        endwhere
        ! Fully attenuated layer detection at subgrid-scale:
-       where ( (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ge. 0.) .and. (x(:,:,k) .ne. undef) ) !DEBUG
+       where ( (x(:,:,k) < S_att_opaq) .and. (x(:,:,k) >= 0.) .and. (x(:,:,k) /= undef) ) !DEBUG
           cldyopaq(:,:,k)=1._wp
        elsewhere
@@ -1358,5 +1358,5 @@
 
        ! Number of usefull sub-column layers:
-       where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
+       where ( (x(:,:,k) > S_att) .and. (x(:,:,k) /= undef) )
           srok(:,:,k)=1._wp
        elsewhere
@@ -1364,5 +1364,5 @@
        endwhere
        ! Number of usefull sub-columns layers for z_opaque 3D fraction:
-       where ( (x(:,:,k) .ge. 0.) .and. (x(:,:,k) .ne. undef) ) !DEBUG
+       where ( (x(:,:,k) >= 0.) .and. (x(:,:,k) /= undef) ) !DEBUG
           srokopaq(:,:,k)=1._wp
        elsewhere
@@ -1397,5 +1397,5 @@
 
      ! Declaring non-opaque cloudy profiles as thin cloud profiles
-	   if ( cldlay(ip,ic,4).gt. 0. .and. cldlay(ip,ic,1) .eq. 0. ) then
+	   if ( cldlay(ip,ic,4)> 0. .and. cldlay(ip,ic,1) == 0. ) then
 	      cldlay(ip,ic,2)  =  1._wp
  	   endif
@@ -1404,5 +1404,5 @@
 
      ! Opaque cloud profiles
-	   if ( cldlay(ip,ic,1) .eq. 1. ) then
+	   if ( cldlay(ip,ic,1) == 1. ) then
 	      zopac = 0._wp
 	      z_top = 0._wp
@@ -1410,5 +1410,5 @@
      ! 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
+	         if ( cldy(ip,ic,Nlevels-k) == 1. .and. zopac == 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
@@ -1416,5 +1416,5 @@
 		    zopac = Nlevels-k+1                        ! z_opaque vertical index on vgrid_z
 		 endif
-	         if ( cldy(ip,ic,Nlevels-k) .eq. 1. ) then
+	         if ( cldy(ip,ic,Nlevels-k) == 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
@@ -1423,5 +1423,5 @@
      ! Summing opaque cloud mean temperatures and altitudes
      ! as defined in Vaillant de Guelis et al. 2017a, AMT
-              if (zopac .ne. 0) then 
+              if (zopac /= 0) then
                  cldtypetemp(ip,1) = cldtypetemp(ip,1) + ( tmp(ip,zopac) + tmp(ip,z_top) )/2.
                  cldtypetemp(ip,3) = cldtypetemp(ip,3) + tmp(ip,zopac)                 ! z_opaque
@@ -1435,5 +1435,5 @@
 
      ! Thin cloud profiles
-	   if ( cldlay(ip,ic,2) .eq. 1. ) then
+	   if ( cldlay(ip,ic,2) == 1. ) then
 	      topcloud = 0._wp
 	      z_top = 0._wp
@@ -1442,9 +1442,9 @@
      ! Declaring thin cloud fraction for 3D variable
      ! From TOA-2-SFC
-                 if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 1. ) then
+                 if ( cldy(ip,ic,k) == 1. .and. topcloud == 1. ) then
                     lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
 		    z_base = k ! bottom cloud layer
                  endif
-	         if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 0. ) then
+	         if ( cldy(ip,ic,k) == 1. .and. topcloud == 0. ) then
                     lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
 		    z_top = k  ! top cloud layer
@@ -1458,5 +1458,5 @@
 	      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
+	         if (  (x(ip,ic,k) > S_att_opaq) .and. (x(ip,ic,k) < 1.0) .and. (x(ip,ic,k) /= undef)  ) then
 		    srmean = srmean + x(ip,ic,k)
 		    srcount = srcount + 1.
@@ -1464,5 +1464,5 @@
 	      enddo
 	      ! If clear sky layers exist below bottom cloud layer
-	      if ( srcount .gt. 0. ) then
+	      if ( srcount > 0. ) then
 	      	 trans2 = srmean/srcount              ! thin cloud transmittance**2
 	      	 tau_app = -(log(trans2))/2.          ! apparent cloud optical depth
@@ -1484,5 +1484,5 @@
 
     ! 3D cloud types fraction (opaque=1 and thin=2 clouds)
-    where ( nsub(:,:) .gt. 0. )
+    where ( nsub(:,:) > 0. )
        lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
        lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
@@ -1492,5 +1492,5 @@
     endwhere
     ! 3D z_opaque fraction (=3)
-    where ( nsubopaq(:,:) .gt. 0. )
+    where ( nsubopaq(:,:) > 0. )
        lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
     elsewhere
@@ -1502,5 +1502,5 @@
     do ip = 1, Npoints
      	do k = 2, Nlevels
-            if ( (lidarcldtype(ip,k,3) .ne. undef) .and. (lidarcldtype(ip,k-1,4) .ne. undef) ) then
+            if ( (lidarcldtype(ip,k,3) /= undef) .and. (lidarcldtype(ip,k-1,4) /= undef) ) then
 	        lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4)
 	    else
@@ -1520,5 +1520,5 @@
 
     ! Mean temperature and altitude
-    where (cldtype(:,1) .gt. 0.)
+    where (cldtype(:,1) > 0.)
        cldtypetemp(:,1) = cldtypetemp(:,1)/cldtype(:,1) ! opaque cloud temp
        cldtypetemp(:,3) = cldtypetemp(:,3)/cldtype(:,1) ! z_opaque
@@ -1534,5 +1534,5 @@
     endwhere
 
-    where (cldtype(:,2) .gt. 0.) ! thin cloud
+    where (cldtype(:,2) > 0.) ! thin cloud
        cldtypetemp(:,2) = cldtypetemp(:,2)/cldtype(:,2)
        cldtypemeanz(:,2) = cldtypemeanz(:,2)/cldtype(:,2)
@@ -1545,5 +1545,5 @@
 
     ! Mean thin cloud emissivity
-    where (count_emis(:) .gt. 0.) ! thin cloud
+    where (count_emis(:) > 0.) ! thin cloud
        cldthinemis(:) = cldthinemis(:)/count_emis(:)
     elsewhere
@@ -1551,5 +1551,5 @@
     endwhere
 
-    where (nsublayer(:,:) .gt. 0.)
+    where (nsublayer(:,:) > 0.)
        cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
     elsewhere
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_interface.F90	(revision 5082)
@@ -275,5 +275,5 @@
         cfg%Lisccp, cfg%Lmisr, cfg%Lmodis, cfg%Lrttov
 
-    if (overlaplmdz.ne.overlap) then
+    if (overlaplmdz/=overlap) then
        print*,'Attention overlaplmdz different de overlap lu dans namelist '
     endif
@@ -282,5 +282,5 @@
    print*,'On passe par using_xios'
  ELSE
-   if (cosp_init_flag .eq. 0) then
+   if (cosp_init_flag == 0) then
 
       ! Initialize the distributional parameters for hydrometeors in radar simulator.
@@ -311,9 +311,9 @@
 
 !!! Ici on modifie les cles logiques pour les outputs selon les champs actives dans les .xml
-  if ((itap.ge.1).and.(first_write))then
+  if ((itap>=1).and.(first_write))then
     IF (using_xios) call read_xiosfieldactive(cfg)
     first_write=.false.
 
-if (cosp_init_flag .eq. 0) then
+if (cosp_init_flag == 0) then
 
     ! Initialize the distributional parameters for hydrometeors in radar simulator.
@@ -384,5 +384,5 @@
 ! 3) Masque terre/mer a partir de la variable fracTerLic
         do ip = 1, Npoints
-          if (fracTerLic(ip).ge.0.5) then
+          if (fracTerLic(ip)>=0.5) then
              land(ip) = 1.
           else
@@ -424,5 +424,5 @@
 
 
-if (cosp_init_flag .eq. 1) then      ! cosp_init_flag = 1
+if (cosp_init_flag == 1) then      ! cosp_init_flag = 1
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -538,5 +538,5 @@
       endif ! debut_cosp
 
-if (cosp_init_flag .eq. 1) then
+if (cosp_init_flag == 1) then
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90	(revision 5082)
@@ -669,6 +669,6 @@
     CHARACTER(LEN=20) :: typeecrit
 
-    ! ug On rÃ©cupÃ¨re le type Ã©crit de la structure:
-    !       Assez moche, Ã|  refaire si meilleure mÃ©thode...
+    ! ug On récupère le type écrit de la structure:
+    !       Assez moche, �|  refaire si meilleure méthode...
     IF (INDEX(var%cosp_typeecrit(iff), "once") > 0) THEN
        typeecrit = 'once'
@@ -730,23 +730,23 @@
 
 ! Axe vertical
-      IF (nvertsave.eq.nvertp(iff)) THEN
+      IF (nvertsave==nvertp(iff)) THEN
           klevs=PARASOL_NREFL
           nam_axvert="sza"
-      ELSE IF (nvertsave.eq.nvertisccp(iff)) THEN
+      ELSE IF (nvertsave==nvertisccp(iff)) THEN
           klevs=7
           nam_axvert="pressure2"
-      ELSE IF (nvertsave.eq.nvertcol(iff)) THEN
+      ELSE IF (nvertsave==nvertcol(iff)) THEN
           klevs=Ncolout
           nam_axvert="column"
-      ELSE IF (nvertsave.eq.nverttemp(iff)) THEN
+      ELSE IF (nvertsave==nverttemp(iff)) THEN
           klevs=LIDAR_NTEMP
           nam_axvert="temp"
-      ELSE IF (nvertsave.eq.nvertmisr(iff)) THEN
+      ELSE IF (nvertsave==nvertmisr(iff)) THEN
           klevs=numMISRHgtBins
           nam_axvert="cth16"
-      ELSE IF (nvertsave.eq.nvertReffIce(iff)) THEN
+      ELSE IF (nvertsave==nvertReffIce(iff)) THEN
           klevs= numMODISReffIceBins
           nam_axvert="ReffIce"
-      ELSE IF (nvertsave.eq.nvertReffLiq(iff)) THEN
+      ELSE IF (nvertsave==nvertReffLiq(iff)) THEN
           klevs= numMODISReffLiqBins
           nam_axvert="ReffLiq"
@@ -765,6 +765,6 @@
       END IF
 
-    ! ug On rÃ©cupÃ¨re le type Ã©crit de la structure:
-    !       Assez moche, Ã|  refaire si meilleure mÃ©thode...
+    ! ug On récupère le type écrit de la structure:
+    !       Assez moche, �|  refaire si meilleure méthode...
     IF (INDEX(var%cosp_typeecrit(iff), "once") > 0) THEN
        typeecrit = 'once'
@@ -827,9 +827,9 @@
     IF (prt_level >= 9) WRITE(lunout,*)'Begin histrwrite2d ',var%name
 
-  ! On regarde si on est dans la phase de dÃ©finition ou d'Ã©criture:
+  ! On regarde si on est dans la phase de définition ou d'écriture:
   IF(.NOT.cosp_varsdefined) THEN
 !$OMP MASTER
       print*,'var, cosp_varsdefined dans cosp_varsdefined ',var%name, cosp_varsdefined
-      !Si phase de dÃ©finition.... on dÃ©finit
+      !Si phase de définition.... on définit
       CALL conf_cospoutputs(var%name,var%cles)
       DO iff=1, 3
@@ -840,5 +840,5 @@
 !$OMP END MASTER
   ELSE
-    !Et sinon on.... Ã©crit
+    !Et sinon on.... écrit
     IF (SIZE(field)/=klon) &
   CALL abort_physic('iophy::histwrite2d_cosp','Field first DIMENSION not equal to klon',1) 
@@ -921,7 +921,7 @@
                nom=var%name
       END IF
-  ! On regarde si on est dans la phase de dÃ©finition ou d'Ã©criture:
+  ! On regarde si on est dans la phase de définition ou d'écriture:
   IF(.NOT.cosp_varsdefined) THEN
-      !Si phase de dÃ©finition.... on dÃ©finit
+      !Si phase de définition.... on définit
 !$OMP MASTER
       CALL conf_cospoutputs(var%name,var%cles)
@@ -933,5 +933,5 @@
 !$OMP END MASTER
   ELSE
-    !Et sinon on.... Ã©crit
+    !Et sinon on.... écrit
     IF (SIZE(field,1)/=klon) &
    CALL abort_physic('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)                                  
@@ -1000,5 +1000,5 @@
 
   IF(cosp_varsdefined) THEN
-    !Et sinon on.... Ã©crit
+    !Et sinon on.... écrit
     IF (SIZE(field,1)/=klon) &
    CALL abort_physic('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)            
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 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90	(revision 5082)
@@ -896,5 +896,5 @@
 !  i = i+1 !si on laisse, 108 au lieu de 107
 
-  if (i.gt.107) then
+  if (i>107) then
      print *, 'COSP_IO: wrong number of output diagnostics'
      print *, i,107
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 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90	(revision 5082)
@@ -98,5 +98,5 @@
     logical :: cmpGases=.true.
 
-    if (Ncolumns .gt. 1) then
+    if (Ncolumns > 1) then
        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ! Generate subcolumns for clouds (SCOPS) and precipitation type (PREC_SCOPS)
@@ -107,5 +107,5 @@
        seed = int(cospstateIN%phalf(:,Nlevels+1))  ! In case of NPoints=1
        ! *NOTE* Chunking will change the seed
-       if (NPoints .gt. 1) seed=int((cospstateIN%phalf(:,Nlevels+1)-minval(cospstateIN%phalf(:,Nlevels+1)))/      &
+       if (NPoints > 1) seed=int((cospstateIN%phalf(:,Nlevels+1)-minval(cospstateIN%phalf(:,Nlevels+1)))/      &
             (maxval(cospstateIN%phalf(:,Nlevels+1))-minval(cospstateIN%phalf(:,Nlevels+1)))*100000) + 1
        call init_rng(rngs, seed)
@@ -145,10 +145,10 @@
           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
-                if (frac_prec(j,i,k) .eq. 1)  prec_ls(j,k) = prec_ls(j,k)+1._wp
-                if (frac_prec(j,i,k) .eq. 2)  prec_cv(j,k) = prec_cv(j,k)+1._wp
-                if (frac_prec(j,i,k) .eq. 3)  prec_cv(j,k) = prec_cv(j,k)+1._wp
-                if (frac_prec(j,i,k) .eq. 3)  prec_ls(j,k) = prec_ls(j,k)+1._wp
+                if (cospIN%frac_out(j,i,k)  == 1)  frac_ls(j,k) = frac_ls(j,k)+1._wp
+                if (cospIN%frac_out(j,i,k)  == 2)  frac_cv(j,k) = frac_cv(j,k)+1._wp
+                if (frac_prec(j,i,k) == 1)  prec_ls(j,k) = prec_ls(j,k)+1._wp
+                if (frac_prec(j,i,k) == 2)  prec_cv(j,k) = prec_cv(j,k)+1._wp
+                if (frac_prec(j,i,k) == 3)  prec_cv(j,k) = prec_cv(j,k)+1._wp
+                if (frac_prec(j,i,k) == 3)  prec_ls(j,k) = prec_ls(j,k)+1._wp
              enddo
              frac_ls(j,k)=frac_ls(j,k)/nColumns
@@ -217,9 +217,9 @@
           do j=1,nPoints
              ! In-cloud mixing ratios.
-             if (frac_ls(j,k) .ne. 0.) then
+             if (frac_ls(j,k) /= 0.) then
                 mr_hydro(j,:,k,I_LSCLIQ) = mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
                 mr_hydro(j,:,k,I_LSCICE) = mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
              endif
-             if (frac_cv(j,k) .ne. 0.) then
+             if (frac_cv(j,k) /= 0.) then
                 mr_hydro(j,:,k,I_CVCLIQ) = mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
                 mr_hydro(j,:,k,I_CVCICE) = mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
@@ -227,20 +227,20 @@
              ! Precipitation
              if (use_precipitation_fluxes) then
-                if (prec_ls(j,k) .ne. 0.) then
+                if (prec_ls(j,k) /= 0.) then
                    fl_lsrain(j,k) = fl_lsrainIN(j,k)/prec_ls(j,k)
                    fl_lssnow(j,k) = fl_lssnowIN(j,k)/prec_ls(j,k)
                    fl_lsgrpl(j,k) = fl_lsgrplIN(j,k)/prec_ls(j,k)
                 endif
-                if (prec_cv(j,k) .ne. 0.) then
+                if (prec_cv(j,k) /= 0.) then
                    fl_ccrain(j,k) = fl_ccrainIN(j,k)/prec_cv(j,k)
                    fl_ccsnow(j,k) = fl_ccsnowIN(j,k)/prec_cv(j,k)
                 endif
              else
-                if (prec_ls(j,k) .ne. 0.) then
+                if (prec_ls(j,k) /= 0.) then
                    mr_hydro(j,:,k,I_LSRAIN) = mr_hydro(j,:,k,I_LSRAIN)/prec_ls(j,k)
                    mr_hydro(j,:,k,I_LSSNOW) = mr_hydro(j,:,k,I_LSSNOW)/prec_ls(j,k)
                    mr_hydro(j,:,k,I_LSGRPL) = mr_hydro(j,:,k,I_LSGRPL)/prec_ls(j,k)
                 endif
-                if (prec_cv(j,k) .ne. 0.) then
+                if (prec_cv(j,k) /= 0.) then
                    mr_hydro(j,:,k,I_CVRAIN) = mr_hydro(j,:,k,I_CVRAIN)/prec_cv(j,k)
                    mr_hydro(j,:,k,I_CVSNOW) = mr_hydro(j,:,k,I_CVSNOW)/prec_cv(j,k)
@@ -361,5 +361,5 @@
        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
+             if (cospIN%rcfg_cloudsat%use_gas_abs == 1 .or. (cospIN%rcfg_cloudsat%use_gas_abs == 2 .and. j == 1)) then
                 g_vol(i,j) = gases(cospstateIN%pfull(i,j), cospstateIN%at(i,j),cospstateIN%qv(i,j),cospIN%rcfg_cloudsat%freq)
              endif
@@ -379,7 +379,7 @@
           
           ! At each model level, what fraction of the precipitation is frozen?
-          where(mr_hydro(:,k,:,I_LSRAIN) .gt. 0 .or. mr_hydro(:,k,:,I_LSSNOW) .gt. 0 .or. &
-                mr_hydro(:,k,:,I_CVRAIN) .gt. 0 .or. mr_hydro(:,k,:,I_CVSNOW) .gt. 0 .or. &
-                mr_hydro(:,k,:,I_LSGRPL) .gt. 0)
+          where(mr_hydro(:,k,:,I_LSRAIN) > 0 .or. mr_hydro(:,k,:,I_LSSNOW) > 0 .or. &
+                mr_hydro(:,k,:,I_CVRAIN) > 0 .or. mr_hydro(:,k,:,I_CVSNOW) > 0 .or. &
+                mr_hydro(:,k,:,I_LSGRPL) > 0)
              fracPrecipIce(:,k,:) = (mr_hydro(:,k,:,I_LSSNOW) + mr_hydro(:,k,:,I_CVSNOW) + &
                   mr_hydro(:,k,:,I_LSGRPL)) / &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mo_rng.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mo_rng.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mo_rng.F90	(revision 5082)
@@ -96,5 +96,5 @@
     !          so we use sizeof(someInt) to determine wheter it is on 32 bit.
     !if ( i2_16*i2_16 .le. huge32 ) then
-    if (digits(testInt) .le. 31) then
+    if (digits(testInt) <= 31) then
     !if (sizeof(testInt) .eq. 4) then
        r=r+1
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/parasol.F90	(revision 5082)
@@ -97,5 +97,5 @@
 
     ! Relative fraction of the opt. thick due to liquid or ice clouds
-    WHERE (tautot_S(1:npoints) .gt. 0.)
+    WHERE (tautot_S(1:npoints) > 0.)
        frac_taucol_liq(1:npoints) = tautot_S_liq(1:npoints) / tautot_S(1:npoints)
        frac_taucol_ice(1:npoints) = tautot_S_ice(1:npoints) / tautot_S(1:npoints)
@@ -118,6 +118,6 @@
     DO it=1,PARASOL_NREFL
        DO ny=1,PARASOL_NTAU-1
-          WHERE (tautot_S(1:npoints) .ge. PARASOL_TAU(ny).and. &
-                 tautot_S(1:npoints) .le. PARASOL_TAU(ny+1))
+          WHERE (tautot_S(1:npoints) >= PARASOL_TAU(ny).and. &
+                 tautot_S(1:npoints) <= PARASOL_TAU(ny+1))
              rlumA_mod(1:npoints,it) = aA(it,ny)*tautot_S(1:npoints) + bA(it,ny)
              rlumB_mod(1:npoints,it) = aB(it,ny)*tautot_S(1:npoints) + bB(it,ny)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90	(revision 5081)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/prec_scops.F90	(revision 5082)
@@ -64,5 +64,5 @@
 
       cv_col = scops_ccfrac*ncol
-      if (cv_col .eq. 0) cv_col=1
+      if (cv_col == 0) cv_col=1
  
       do ilev=1,nlev
@@ -81,15 +81,15 @@
         flag_cv=0
         do ilev=1,nlev
-          if (frac_out(j,ibox,ilev) .eq. 1) then 
+          if (frac_out(j,ibox,ilev) == 1) then
             flag_ls=1
           endif
-          if (frac_out(j,ibox,ilev) .eq. 2) then 
+          if (frac_out(j,ibox,ilev) == 2) then
             flag_cv=1
           endif
         enddo !loop over nlev
-        if (flag_ls .eq. 1) then
+        if (flag_ls == 1) then
            frac_out_ls(j,ibox)=1
         endif
-        if (flag_cv .eq. 1) then
+        if (flag_cv == 1) then
            frac_out_cv(j,ibox)=1
         endif
@@ -102,14 +102,14 @@
         flag_cv=0
     
-        if (ls_p_rate(j,1) .gt. 0.) then 
+        if (ls_p_rate(j,1) > 0.) then
             do ibox=1,ncol ! possibility ONE
-                if (frac_out(j,ibox,1) .eq. 1) then 
+                if (frac_out(j,ibox,1) == 1) then
                     prec_frac(j,ibox,1) = 1
                     flag_ls=1
                 endif
             enddo ! loop over ncol
-            if (flag_ls .eq. 0) then ! possibility THREE
+            if (flag_ls == 0) then ! possibility THREE
                 do ibox=1,ncol
-                    if (frac_out(j,ibox,2) .eq. 1) then 
+                    if (frac_out(j,ibox,2) == 1) then
                         prec_frac(j,ibox,1) = 1
                         flag_ls=1
@@ -117,7 +117,7 @@
                 enddo ! loop over ncol
             endif
-        if (flag_ls .eq. 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_ls(j,ibox) .eq. 1) then 
+        if (flag_ls == 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_ls(j,ibox) == 1) then
             prec_frac(j,ibox,1) = 1
             flag_ls=1
@@ -125,5 +125,5 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls .eq. 0) then ! possibility Five
+        if (flag_ls == 0) then ! possibility Five
         do ibox=1,ncol
     !     prec_frac(j,1:ncol,1) = 1
@@ -134,8 +134,8 @@
        ! There is large scale precipitation
      
-        if (cv_p_rate(j,1) .gt. 0.) then 
+        if (cv_p_rate(j,1) > 0.) then
          do ibox=1,ncol ! possibility ONE
-          if (frac_out(j,ibox,1) .eq. 2) then 
-           if (prec_frac(j,ibox,1) .eq. 0) then
+          if (frac_out(j,ibox,1) == 2) then
+           if (prec_frac(j,ibox,1) == 0) then
         prec_frac(j,ibox,1) = 2
        else
@@ -145,8 +145,8 @@
       endif
         enddo ! loop over ncol
-        if (flag_cv .eq. 0) then ! possibility THREE
-        do ibox=1,ncol
-        if (frac_out(j,ibox,2) .eq. 2) then 
-                if (prec_frac(j,ibox,1) .eq. 0) then
+        if (flag_cv == 0) then ! possibility THREE
+        do ibox=1,ncol
+        if (frac_out(j,ibox,2) == 2) then
+                if (prec_frac(j,ibox,1) == 0) then
             prec_frac(j,ibox,1) = 2
             else
@@ -157,8 +157,8 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv .eq. 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_cv(j,ibox) .eq. 1) then 
-                if (prec_frac(j,ibox,1) .eq. 0) then
+        if (flag_cv == 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_cv(j,ibox) == 1) then
+                if (prec_frac(j,ibox,1) == 0) then
             prec_frac(j,ibox,1) = 2
             else
@@ -169,7 +169,7 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv .eq. 0) then  ! possibility Five
+        if (flag_cv == 0) then  ! possibility Five
         do ibox=1,cv_col
-                if (prec_frac(j,ibox,1) .eq. 0) then
+                if (prec_frac(j,ibox,1) == 0) then
             prec_frac(j,ibox,1) = 2
             else
@@ -192,15 +192,15 @@
         flag_cv=0
     
-        if (ls_p_rate(j,ilev) .gt. 0.) then 
+        if (ls_p_rate(j,ilev) > 0.) then
          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 
+          if ((frac_out(j,ibox,ilev) == 1) .or. ((prec_frac(j,ibox,ilev-1) == 1)     &
+            .or. (prec_frac(j,ibox,ilev-1) == 3))) then
            prec_frac(j,ibox,ilev) = 1
            flag_ls=1
           endif
         enddo ! loop over ncol
-        if ((flag_ls .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
-        do ibox=1,ncol
-        if (frac_out(j,ibox,ilev+1) .eq. 1) then 
+        if ((flag_ls == 0) .and. (ilev < nlev)) then ! possibility THREE
+        do ibox=1,ncol
+        if (frac_out(j,ibox,ilev+1) == 1) then
             prec_frac(j,ibox,ilev) = 1
             flag_ls=1
@@ -208,7 +208,7 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls .eq. 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_ls(j,ibox) .eq. 1) then 
+        if (flag_ls == 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_ls(j,ibox) == 1) then
             prec_frac(j,ibox,ilev) = 1
             flag_ls=1
@@ -216,5 +216,5 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls .eq. 0) then ! possibility Five
+        if (flag_ls == 0) then ! possibility Five
         do ibox=1,ncol
 !     prec_frac(j,1:ncol,ilev) = 1
@@ -224,9 +224,9 @@
       endif ! There is large scale precipitation
     
-        if (cv_p_rate(j,ilev) .gt. 0.) then 
+        if (cv_p_rate(j,ilev) > 0.) then
          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 
-            if (prec_frac(j,ibox,ilev) .eq. 0) then
+          if ((frac_out(j,ibox,ilev) == 2) .or. ((prec_frac(j,ibox,ilev-1) == 2)     &
+            .or. (prec_frac(j,ibox,ilev-1) == 3))) then
+            if (prec_frac(j,ibox,ilev) == 0) then
          prec_frac(j,ibox,ilev) = 2
         else
@@ -236,8 +236,8 @@
         endif
        enddo ! loop over ncol
-        if ((flag_cv .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
-        do ibox=1,ncol
-        if (frac_out(j,ibox,ilev+1) .eq. 2) then 
-                if (prec_frac(j,ibox,ilev) .eq. 0) then
+        if ((flag_cv == 0) .and. (ilev < nlev)) then ! possibility THREE
+        do ibox=1,ncol
+        if (frac_out(j,ibox,ilev+1) == 2) then
+                if (prec_frac(j,ibox,ilev) == 0) then
             prec_frac(j,ibox,ilev) = 2
             else
@@ -248,8 +248,8 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv .eq. 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_cv(j,ibox) .eq. 1) then 
-                if (prec_frac(j,ibox,ilev) .eq. 0) then
+        if (flag_cv == 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_cv(j,ibox) == 1) then
+                if (prec_frac(j,ibox,ilev) == 0) then
             prec_frac(j,ibox,ilev) = 2
             else
@@ -260,7 +260,7 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv .eq. 0) then  ! possibility Five 
+        if (flag_cv == 0) then  ! possibility Five
         do ibox=1,cv_col
-                if (prec_frac(j,ibox,ilev) .eq. 0) then
+                if (prec_frac(j,ibox,ilev) == 0) then
             prec_frac(j,ibox,ilev) = 2
             else
