Changeset 5095 for LMDZ6/branches/Amaury_dev/libf/phylmd/cosp

Timestamp:

Jul 22, 2024, 9:46:57 AM (12 months ago)

Author:

abarral

Message:

Revert cosp*/ from the trunk, as it's external code
Add missing bits from FCM2 source

Location:

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp

Files:

• MISR_simulator.F (modified) (17 diffs)
• array_lib.F90 (modified) (2 diffs)
• atmos_lib.F90 (modified) (1 diff)
• calc_Re.F90 (modified) (5 diffs)
• cosp_output_write_mod.F90 (modified) (9 diffs)
• dsd.F90 (modified) (3 diffs)
• format_input.F90 (modified) (2 diffs)
• gases.F90 (modified) (1 diff)
• icarus.F (modified) (51 diffs)
• lidar_simulator.F90 (modified) (20 diffs)
• m_mrgrnk.F90 (modified) (2 diffs)
• math_lib.F90 (modified) (6 diffs)
• mod_cosp.F90 (modified) (4 diffs)
• mod_cosp_modis_simulator.F90 (modified) (3 diffs)
• mod_cosp_radar.F90 (modified) (3 diffs)
• mod_cosp_types.F90 (modified) (1 diff)
• mod_llnl_stats.F90 (modified) (2 diffs)
• mod_lmd_ipsl_stats.F90 (modified) (51 diffs)
• mod_modis_sim.F90 (modified) (9 diffs)
• optics_lib.F90 (modified) (10 diffs)
• pf_to_mr.F (modified) (2 diffs)
• phys_cosp.F90 (modified) (3 diffs)
• prec_scops.F (modified) (16 diffs)
• predict_mom07.F90 (modified) (2 diffs)
• radar_simulator.F90 (modified) (8 diffs)
• radar_simulator_init.F90 (modified) (1 diff)
• radar_simulator_types.F90 (modified) (3 diffs)
• scops.F (modified) (13 diffs)
• zeff.F90 (modified) (3 diffs)

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/MISR_simulator.F

@@ -133 +133 @@
             
         ! define location of "layer top"
-        if(ilev==1 .or. ilev==nlev) then
+        if(ilev.eq.1 .or. ilev.eq.nlev) then
             ztest=zfull(j,ilev)
         else
@@ -144 +144 @@
         do loop=2,n_MISR_CTH
         
-            if ( ztest >
+            if ( ztest .gt.
      &                1000*MISR_CTH_boundaries(loop+1) ) then
         
@@ -173 +173 @@
              dtau=0
              
-             if (frac_out(j,ibox,ilev)==1) then
+             if (frac_out(j,ibox,ilev).eq.1) then
                         dtau = dtau_s(j,ilev)
                  endif
                  
-                 if (frac_out(j,ibox,ilev)==2) then
+                 if (frac_out(j,ibox,ilev).eq.2) then
                         dtau = dtau_c(j,ilev)
                  end if 
@@ -186 +186 @@
         ! NOW for MISR ..
         ! if there a cloud ... start the counter ... store this height
-        if(thres_crossed_MISR == 0 .and. dtau > 0.) then
+        if(thres_crossed_MISR .eq. 0 .and. dtau .gt. 0.) then
         
             ! first encountered a "cloud"
@@ -193 +193 @@
         endif   
                 
-        if( thres_crossed_MISR < 99 .and.
-     &              thres_crossed_MISR > 0 ) then
+        if( thres_crossed_MISR .lt. 99 .and.
+     &              thres_crossed_MISR .gt. 0 ) then
      
-                if( dtau == 0.) then
+                if( dtau .eq. 0.) then
         
                     ! we have come to the end of the current cloud
@@ -212 +212 @@
             ! then MISR will like see a top below the top of the current 
             ! layer
-            if( dtau>0 .and. (cloud_dtau-dtau) < 1) then
-            
-                if(dtau < 1 .or. ilev==1 .or. ilev==nlev) then
+            if( dtau.gt.0 .and. (cloud_dtau-dtau) .lt. 1) then
+            
+                if(dtau .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then
 
                     ! MISR will likely penetrate to some point
@@ -233 +233 @@
         
             ! check for a distinctive water layer
-            if(dtau > 1 .and. at(j,ilev)>273 ) then
+            if(dtau .gt. 1 .and. at(j,ilev).gt.273 ) then
      
                     ! must be a water cloud ... 
@@ -242 +242 @@
             ! if the total column optical depth is "large" than
             ! MISR can't seen anything else ... set current point as CTH level
-            if(tau(j,ibox) > 5) then
+            if(tau(j,ibox) .gt. 5) then 
 
                 thres_crossed_MISR=99           
@@ -254 +254 @@
         ! check to see if there was a cloud for which we didn't 
         ! set a MISR cloud top boundary
-        if( thres_crossed_MISR == 1) then
+        if( thres_crossed_MISR .eq. 1) then
     
         ! if the cloud has a total optical depth of greater
@@ -260 +260 @@
         ! with a height near the true cloud top
         ! otherwise there should be no CTH
-        if( tau(j,ibox) > 0.5) then
+        if( tau(j,ibox) .gt. 0.5) then
 
             ! keep MISR detected CTH
             
-        elseif(tau(j,ibox) > 0.2) then
+        elseif(tau(j,ibox) .gt. 0.2) then
 
             ! MISR may detect but wont likley have a good height
@@ -294 +294 @@
     !   This setup assumes the columns represent a about a 1 to 4 km scale
     !   it will need to be modified significantly, otherwise
-        if(ncol==1) then
+        if(ncol.eq.1) then
     
        ! adjust based on neightboring points ... i.e. only 2D grid was input
            do j=2,npoints-1
             
-            if(box_MISR_ztop(j-1,1)>0 .and.
-     &             box_MISR_ztop(j+1,1)>0       ) then
+            if(box_MISR_ztop(j-1,1).gt.0 .and. 
+     &             box_MISR_ztop(j+1,1).gt.0       ) then
 
                 if( abs( box_MISR_ztop(j-1,1) -  
-     &                   box_MISR_ztop(j+1,1) ) < 500
+     &                   box_MISR_ztop(j+1,1) ) .lt. 500 
      &              .and.
-     &                   box_MISR_ztop(j,1) <
+     &                   box_MISR_ztop(j,1) .lt. 
      &                   box_MISR_ztop(j+1,1)     ) then
             
@@ -319 +319 @@
          do ibox=2,ncol-1
             
-            if(box_MISR_ztop(1,ibox-1)>0 .and.
-     &             box_MISR_ztop(1,ibox+1)>0        ) then
+            if(box_MISR_ztop(1,ibox-1).gt.0 .and. 
+     &             box_MISR_ztop(1,ibox+1).gt.0        ) then
 
                 if( abs( box_MISR_ztop(1,ibox-1) -  
-     &                   box_MISR_ztop(1,ibox+1) ) < 500
+     &                   box_MISR_ztop(1,ibox+1) ) .lt. 500 
      &              .and.
-     &                   box_MISR_ztop(1,ibox) <
+     &                   box_MISR_ztop(1,ibox) .lt. 
      &                   box_MISR_ztop(1,ibox+1)     ) then
             
@@ -357 +357 @@
          do ibox=1,ncol
 
-            if (tau(j,ibox) > (tauchk)) then
+            if (tau(j,ibox) .gt. (tauchk)) then
                box_cloudy(j,ibox)=.true.
             endif
@@ -366 +366 @@
     
           !determine optical depth category
-              if (tau(j,ibox) < isccp_taumin) then
+              if (tau(j,ibox) .lt. isccp_taumin) then
                   itau=1
-              else if (tau(j,ibox) >= isccp_taumin
-     &          .and. tau(j,ibox) < 1.3) then
+              else if (tau(j,ibox) .ge. isccp_taumin                                    
+     &          .and. tau(j,ibox) .lt. 1.3) then
                   itau=2
-              else if (tau(j,ibox) >= 1.3
-     &          .and. tau(j,ibox) < 3.6) then
+              else if (tau(j,ibox) .ge. 1.3 
+     &          .and. tau(j,ibox) .lt. 3.6) then
                   itau=3
-              else if (tau(j,ibox) >= 3.6
-     &          .and. tau(j,ibox) < 9.4) then
+              else if (tau(j,ibox) .ge. 3.6 
+     &          .and. tau(j,ibox) .lt. 9.4) then
                   itau=4
-              else if (tau(j,ibox) >= 9.4
-     &          .and. tau(j,ibox) < 23.) then
+              else if (tau(j,ibox) .ge. 9.4 
+     &          .and. tau(j,ibox) .lt. 23.) then
                   itau=5
-              else if (tau(j,ibox) >= 23.
-     &          .and. tau(j,ibox) < 60.) then
+              else if (tau(j,ibox) .ge. 23. 
+     &          .and. tau(j,ibox) .lt. 60.) then
                   itau=6
-              else if (tau(j,ibox) >= 60.) then
+              else if (tau(j,ibox) .ge. 60.) then
                   itau=7
               endif
@@ -390 +390 @@
 
        ! update MISR histograms and summary metrics - roj 5/2005
-       if (sunlit(j)==1) then
+       if (sunlit(j).eq.1) then 
                      
               !if cloudy added by roj 5/2005
-          if( box_MISR_ztop(j,ibox)==0) then
+          if( box_MISR_ztop(j,ibox).eq.0) then
           
             ! no cloud detected
             iMISR_ztop=0
 
-          elseif( box_MISR_ztop(j,ibox)==-1) then
+          elseif( box_MISR_ztop(j,ibox).eq.-1) then
 
             ! cloud can be detected but too thin to get CTH
@@ -416 +416 @@
             do loop=2,n_MISR_CTH
         
-                if ( box_MISR_ztop(j,ibox) >
+                if ( box_MISR_ztop(j,ibox) .gt.
      &                1000*MISR_CTH_boundaries(loop+1) ) then
         
@@ -466 +466 @@
        enddo ! ibox - loop over subcolumns          
       
-       if( MISR_cldarea(j) > 0.) then
+       if( MISR_cldarea(j) .gt. 0.) then
         MISR_mean_ztop(j)= MISR_mean_ztop(j) / MISR_cldarea(j)   ! roj 5/2006
        endif

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/array_lib.F90

@@ -45 +45 @@
  
 ! ----- INPUTS -----
-  real(kind=8), dimension(:), intent(in) :: list
-  real(kind=8), intent(in) :: val
+  real*8, dimension(:), intent(in) :: list
+  real*8, intent(in) :: val  
   integer, intent(in), optional :: sort
   
 ! ----- OUTPUTS -----
-  integer(kind=4) :: infind
-  real(kind=8), intent(out), optional :: dist
+  integer*4 :: infind
+  real*8, intent(out), optional :: dist
 
 ! ----- INTERNAL -----
-  real(kind=8), dimension(size(list)) :: lists
-  integer(kind=4) :: nlist, result, tmp(1), sort_list
-  integer(kind=4), dimension(size(list)) :: mask, idx
+  real*8, dimension(size(list)) :: lists
+  integer*4 :: nlist, result, tmp(1), sort_list
+  integer*4, dimension(size(list)) :: mask, idx
 
   if (present(sort)) then
@@ -121 +121 @@
 
 ! ----- INPUTS -----
-  real(kind=8), dimension(:), intent(in) :: yarr, xarr, xxarr
-  real(kind=8), intent(in) :: tol
+  real*8, dimension(:), intent(in) :: yarr, xarr, xxarr
+  real*8, intent(in) :: tol
 
 ! ----- OUTPUTS -----
-  real(kind=8), dimension(size(xxarr)), intent(out) :: yyarr
+  real*8, dimension(size(xxarr)), intent(out) :: yyarr
 
 ! ----- INTERNAL -----
-  real(kind=8), dimension(size(xarr)) :: ysort, xsort
-  integer(kind=4), dimension(size(xarr)) :: ist
-  integer(kind=4) :: nx, nxx, i, iloc
-  real(kind=8) :: d, m
+  real*8, dimension(size(xarr)) :: ysort, xsort
+  integer*4, dimension(size(xarr)) :: ist
+  integer*4 :: nx, nxx, i, iloc
+  real*8 :: d, m
 
   nx = size(xarr)

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/atmos_lib.F90

@@ -43 +43 @@
 
 ! ----- OUTPUTS -----
-  real(kind=8), intent(out), dimension(ndat) :: &
+  real*8, intent(out), dimension(ndat) :: &
   hgt, &                        ! height (m)
   prs, &                        ! pressure (hPa)

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/calc_Re.F90

@@ -40 +40 @@
 ! ----- INPUTS -----  
   
-  real(kind=8), intent(in) :: Q,Np,rho_a
+  real*8, intent(in) :: Q,Np,rho_a
  
   integer, intent(in):: dtype
-  real(kind=8), intent(in) :: dmin,dmax,rho_c,p1,p2,p3
-    
-  real(kind=8), intent(inout) :: apm,bpm
+  real*8, intent(in) :: dmin,dmax,rho_c,p1,p2,p3
+    
+  real*8, intent(inout) :: apm,bpm  
     
 ! ----- OUTPUTS -----
 
-  real(kind=8), intent(out) :: Re
+  real*8, intent(out) :: Re
   
 ! ----- INTERNAL -----
  
   integer :: local_dtype
-  real(kind=8)  :: local_p3,local_Np
-
-  real(kind=8) :: pi, &
+  real*8  :: local_p3,local_Np
+
+  real*8 :: pi, &
   N0,D0,vu,dm,ld, &         ! gamma, exponential variables
   rg,log_sigma_g
   
-  real(kind=8) :: tmp1,tmp2
+  real*8 :: tmp1,tmp2
 
   pi = acos(-1.0)
@@ -72 +72 @@
   ! Exponential is same as modified gamma with vu =1
   ! if Np is specified then we will just treat as modified gamma
-  if(dtype==2 .and. Np>0) then
+  if(dtype.eq.2 .and. Np>0) then
     local_dtype=1;
     local_p3=1;
@@ -119 +119 @@
     
 
-    if( Np==0 .and. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default
+    if( Np.eq.0 .and. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default
       
         dm = p2             ! by definition, should have units of microns
@@ -126 +126 @@
     else   ! use value of Np
         
-        if(Np==0) then
+        if(Np.eq.0) then
         
             if( abs(p1+1) > 1E-8 ) then  !   use default number concentration 
@@ -233 +233 @@
      
     ! get rg ... 
-    if( Np==0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
+    if( Np.eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
     
             rg = p2

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_output_write_mod.F90

@@ -185 +185 @@
     do k=1,PARASOL_NREFL
      do ip=1, Npoints
-      if (stlidar%cldlayer(ip,4)>0.01.and.stlidar%parasolrefl(ip,k)/=missing_val) then
+      if (stlidar%cldlayer(ip,4).gt.0.01.and.stlidar%parasolrefl(ip,k).ne.missing_val) then
         parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ &
                              stlidar%cldlayer(ip,4)
@@ -456 +456 @@
     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'
@@ -523 +523 @@
 
 ! Axe vertical
-      IF (nvertsave==nvertp(iff)) THEN
+      IF (nvertsave.eq.nvertp(iff)) THEN
           klevs=PARASOL_NREFL
           nam_axvert="sza"
-      ELSE IF (nvertsave==nvertisccp(iff)) THEN
+      ELSE IF (nvertsave.eq.nvertisccp(iff)) THEN
           klevs=7
           nam_axvert="pressure2"
-      ELSE IF (nvertsave==nvertcol(iff)) THEN
+      ELSE IF (nvertsave.eq.nvertcol(iff)) THEN
           klevs=Ncolout
           nam_axvert="column"
-      ELSE IF (nvertsave==nverttemp(iff)) THEN
+      ELSE IF (nvertsave.eq.nverttemp(iff)) THEN
           klevs=LIDAR_NTEMP
           nam_axvert="temp"
-      ELSE IF (nvertsave==nvertmisr(iff)) THEN
+      ELSE IF (nvertsave.eq.nvertmisr(iff)) THEN
           klevs=MISR_N_CTH
           nam_axvert="cth16"
-      ELSE IF (nvertsave==nvertReffIce(iff)) THEN
+      ELSE IF (nvertsave.eq.nvertReffIce(iff)) THEN
           klevs= numMODISReffIceBins
           nam_axvert="ReffIce"
-      ELSE IF (nvertsave==nvertReffLiq(iff)) THEN
+      ELSE IF (nvertsave.eq.nvertReffLiq(iff)) THEN
           klevs= numMODISReffLiqBins
           nam_axvert="ReffLiq"
@@ -558 +558 @@
       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'
@@ -628 +628 @@
     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
-      !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
@@ -640 +640 @@
 !$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) 
@@ -725 +725 @@
                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)
@@ -737 +737 @@
 !$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)                                  
@@ -809 +809 @@
 
   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)            

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/dsd.F90

@@ -60 +60 @@
   integer, intent(in) :: nsizes
   integer, intent(in) :: dtype
-  real(kind=8), intent(in)  :: Q,Re_,Np,D(nsizes)
-  real(kind=8), intent(in)  :: rho_a,tk,dmin,dmax,rho_c,p1,p2,p3
-    
-  real(kind=8), intent(inout) :: apm,bpm
+  real*8, intent(in)  :: Q,Re_,Np,D(nsizes)
+  real*8, intent(in)  :: rho_a,tk,dmin,dmax,rho_c,p1,p2,p3
+    
+  real*8, intent(inout) :: apm,bpm  
   
 ! ----- OUTPUTS -----
 
-  real(kind=8), intent(out) :: N(nsizes)
+  real*8, intent(out) :: N(nsizes)
   
 ! ----- INTERNAL -----
  
-   real(kind=8) :: fc(nsizes)
-
-  real(kind=8) :: &
+   real*8 :: fc(nsizes)
+
+  real*8 :: &
   N0,D0,vu,local_np,dm,ld, &            ! gamma, exponential variables
   dmin_mm,dmax_mm,ahp,bhp, &        ! power law variables
@@ -79 +79 @@
   rho_e                 ! particle density (kg m^-3)
   
-  real(kind=8) :: tmp1, tmp2
-  real(kind=8) :: pi,rc,tc
-  real(kind=8) :: Re
+  real*8 :: tmp1, tmp2
+  real*8 :: pi,rc,tc
+  real*8 :: Re
 
   integer k,lidx,uidx
@@ -352 +352 @@
       log_sigma_g = p3
       tmp2 = (bpm*log_sigma_g)**2.
-      if(Re<=0) then
+      if(Re.le.0) then 
         rg = p2
       else

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/format_input.F90

@@ -40 +40 @@
 
 ! ----- INPUTS -----
-  real(kind=8), dimension(:,:), intent(in) :: &
+  real*8, dimension(:,:), intent(in) :: &
     hgt_matrix,env_hgt_matrix,env_t_matrix,env_p_matrix,env_rh_matrix
 
 ! ----- OUTPUTS -----
-  real(kind=8), dimension(:,:), intent(out) :: &
+  real*8, dimension(:,:), intent(out) :: &
     t_matrix,p_matrix,rh_matrix
 
@@ -97 +97 @@
 
 ! ----- OUTPUTS -----
-  real(kind=8), dimension(:,:), intent(inout) :: &
+  real*8, dimension(:,:), intent(inout) :: &
     hgt_matrix,p_matrix,t_matrix,rh_matrix
-  real(kind=8), dimension(:,:,:), intent(inout) :: &
+  real*8, dimension(:,:,:), intent(inout) :: &
     hm_matrix
   logical, intent(out) :: hgt_reversed

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/gases.F90

@@ -27 +27 @@
   nbands_o2 = 48 ,&
   nbands_h2o = 30
-  real(kind=8), intent(in) :: PRES_mb, T, RH, f
-  real(kind=8) :: gases, th, e, p, sumo, gm0, a0, ap, term1, term2, term3, &
+  real*8, intent(in) :: PRES_mb, T, RH, f
+  real*8 :: gases, th, e, p, sumo, gm0, a0, ap, term1, term2, term3, &
             bf, be, term4, npp
-  real(kind=8), dimension(nbands_o2) :: v0, a1, a2, a3, a4, a5, a6
-  real(kind=8), dimension(nbands_h2o) :: v1, b1, b2, b3
-  real(kind=8) :: e_th,one_th,pth3,eth35,aux1,aux2,aux3,aux4
-  real(kind=8) :: gm,delt,x,y,gm2
-  real(kind=8) :: fpp_o2,fpp_h2o,s_o2,s_h2o
+  real*8, dimension(nbands_o2) :: v0, a1, a2, a3, a4, a5, a6
+  real*8, dimension(nbands_h2o) :: v1, b1, b2, b3
+  real*8 :: e_th,one_th,pth3,eth35,aux1,aux2,aux3,aux4
+  real*8 :: gm,delt,x,y,gm2
+  real*8 :: fpp_o2,fpp_h2o,s_o2,s_h2o
   integer :: i
   

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/icarus.F

@@ -293 +293 @@
       ncolprint=0
 
-      if ( debug/=0 ) then
+      if ( debug.ne.0 ) then
           j=1
           write(6,'(a10)') 'j='
@@ -347 +347 @@
 !     ---------------------------------------------------!
 
-      if (ncolprint/=0) then
+      if (ncolprint.ne.0) then
       do j=1,npoints,1000
         write(6,'(a10)') 'j='
@@ -354 +354 @@
       endif
 
-      if (top_height == 1 .or. top_height == 3) then
+      if (top_height .eq. 1 .or. top_height .eq. 3) then 
 
       do j=1,npoints 
@@ -364 +364 @@
       enddo 
 
-      do ilev=1,nlev
+      do 12 ilev=1,nlev
         do j=1,npoints 
-         if (pfull(j,ilev) < 40000. .and.
-     &          pfull(j,ilev) >  5000. .and.
-     &          at(j,ilev) < attropmin(j)) then
+         if (pfull(j,ilev) .lt. 40000. .and.
+     &          pfull(j,ilev) .gt.  5000. .and.
+     &          at(j,ilev) .lt. attropmin(j)) then
                 ptrop(j) = pfull(j,ilev)
                 attropmin(j) = at(j,ilev)
@@ -375 +375 @@
            end if
         enddo
-      END DO
-
-      do ilev=1,nlev
+12    continue
+
+      do 13 ilev=1,nlev
         do j=1,npoints 
-           if (at(j,ilev) > atmax(j) .and.
-     &              ilev  >= itrop(j)) atmax(j)=at(j,ilev)
-        enddo
-      END DO
+           if (at(j,ilev) .gt. atmax(j) .and.
+     &              ilev  .ge. itrop(j)) atmax(j)=at(j,ilev)
+        enddo
+13    continue
 
       end if
 
 
-      if (top_height == 1 .or. top_height == 3) then
+      if (top_height .eq. 1 .or. top_height .eq. 3) then
           do j=1,npoints
               meantb(j) = 0.
               meantbclr(j) = 0. 
-          END DO
+          end do
       else
           do j=1,npoints
               meantb(j) = output_missing_value
               meantbclr(j) = output_missing_value
-          END DO
+          end do
       end if
       
@@ -408 +408 @@
           rangevec(j)=0
 
-          if (cc(j,ilev) < 0. .or. cc(j,ilev) > 1.) then
+          if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
 !           error = cloud fraction less than zero
 !           error = cloud fraction greater than 1
@@ -414 +414 @@
           endif
 
-          if (conv(j,ilev) < 0. .or. conv(j,ilev) > 1.) then
+          if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
 !           ' error = convective cloud fraction less than zero'
 !           ' error = convective cloud fraction greater than 1'
@@ -420 +420 @@
           endif
 
-          if (dtau_s(j,ilev) < 0.) then
+          if (dtau_s(j,ilev) .lt. 0.) then
 !           ' error = stratiform cloud opt. depth less than zero'
             rangevec(j)=rangevec(j)+4
           endif
 
-          if (dtau_c(j,ilev) < 0.) then
+          if (dtau_c(j,ilev) .lt. 0.) then
 !           ' error = convective cloud opt. depth less than zero'
             rangevec(j)=rangevec(j)+8
           endif
 
-          if (dem_s(j,ilev) < 0. .or. dem_s(j,ilev) > 1.) then
+          if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
 !             ' error = stratiform cloud emissivity less than zero'
 !             ' error = stratiform cloud emissivity greater than 1'
@@ -436 +436 @@
           endif
 
-          if (dem_c(j,ilev) < 0. .or. dem_c(j,ilev) > 1.) then
+          if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
 !             ' error = convective cloud emissivity less than zero'
 !             ' error = convective cloud emissivity greater than 1'
@@ -448 +448 @@
         enddo
 
-        if (rangeerror/=0) then
+        if (rangeerror.ne.0) then 
               write (6,*) 'Input variable out of range'
               write (6,*) 'rangevec:'
@@ -466 +466 @@
  
       !initialize tau and albedocld to zero
-      do ibox=1,ncol
+      do 15 ibox=1,ncol
         do j=1,npoints 
             tau(j,ibox)=0.
@@ -474 +474 @@
           box_cloudy(j,ibox)=.false.
         enddo
-      END DO
+15    continue
 
       !compute total cloud optical depth for each column     
@@ -481 +481 @@
             do ibox=1,ncol
               do j=1,npoints 
-                 if (frac_out(j,ibox,ilev)==1) then
+                 if (frac_out(j,ibox,ilev).eq.1) then
                         tau(j,ibox)=tau(j,ibox)
      &                     + dtau_s(j,ilev)
                  endif
-                 if (frac_out(j,ibox,ilev)==2) then
+                 if (frac_out(j,ibox,ilev).eq.2) then
                         tau(j,ibox)=tau(j,ibox)
      &                     + dtau_c(j,ilev)
@@ -492 +492 @@
             enddo ! ibox
       enddo ! ilev
-          if (ncolprint/=0) then
+          if (ncolprint.ne.0) then
 
               do j=1,npoints ,1000
@@ -521 +521 @@
 !             sky versions of these quantities.
 
-      if (top_height == 1 .or. top_height == 3) then
+      if (top_height .eq. 1 .or. top_height .eq. 3) then
 
 
@@ -539 +539 @@
         pstd = 1.013250E+06
         t0 = 296.
-        if (ncolprint /= 0)
+        if (ncolprint .ne. 0) 
      &         write(6,*)  'ilev   pw (kg/m2)   tauwv(j)      dem_wv'
-        do ilev=1,nlev
+        do 125 ilev=1,nlev
           do j=1,npoints 
                !press and dpress are dyne/cm2 = Pascals *10
@@ -559 +559 @@
                dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
           enddo
-               if (ncolprint /= 0) then
+               if (ncolprint .ne. 0) then
                do j=1,npoints ,1000
                write(6,'(a10)') 'j='
@@ -568 +568 @@
                enddo
              endif
-      END DO
+125     continue
 
         !initialize variables
@@ -598 +598 @@
 
           enddo   
-            if (ncolprint/=0) then
+            if (ncolprint.ne.0) then
              do j=1,npoints ,1000
               write(6,'(a10)') 'j='
@@ -627 +627 @@
         enddo
 
-        if (ncolprint/=0) then
+        if (ncolprint.ne.0) then
         do j=1,npoints ,1000
           write(6,'(a10)') 'j='
@@ -649 +649 @@
 
 
-        if (ncolprint/=0) then
+        if (ncolprint.ne.0) then
 
         do j=1,npoints ,1000
@@ -683 +683 @@
 
               ! emissivity for point in this layer
-                if (frac_out(j,ibox,ilev)==1) then
+                if (frac_out(j,ibox,ilev).eq.1) then
                 dem(j,ibox)= 1. - 
      &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_s(j,ilev)) )
-                else if (frac_out(j,ibox,ilev)==2) then
+                else if (frac_out(j,ibox,ilev).eq.2) then
                 dem(j,ibox)= 1. - 
      &             ( (1. - dem_wv(j,ilev)) * (1. -  dem_c(j,ilev)) )
@@ -710 +710 @@
             enddo ! ibox
 
-            if (ncolprint/=0) then
+            if (ncolprint.ne.0) then
               do j=1,npoints,1000
               write (6,'(a)') 'ilev:'
@@ -740 +740 @@
             bb(j)=1/( exp(1307.27/skt(j)) - 1. )
             !bb(j)=5.67e-8*skt(j)**4
-          END DO
+          end do
 
         do ibox=1,ncol
@@ -751 +751 @@
      &         * trans_layers_above(j,ibox) 
             
-          END DO
-        END DO
+          end do
+        end do
 
         !calculate mean infrared brightness temperature
@@ -758 +758 @@
           do j=1,npoints 
             meantb(j) = meantb(j)+1307.27/(log(1.+(1./fluxtop(j,ibox))))
-          END DO
-        END DO
+          end do
+        end do
           do j=1, npoints
             meantb(j) = meantb(j) / real(ncol)
-          END DO
-
-        if (ncolprint/=0) then
+          end do        
+
+        if (ncolprint.ne.0) then
 
           do j=1,npoints ,1000
@@ -784 +784 @@
           write (6,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
       
-          END DO
+          end do
       endif
     
@@ -819 +819 @@
           enddo 
 
-          if (top_height == 1) then
+          if (top_height .eq. 1) then
             do j=1,npoints  
-              if (transmax(j) > 0.001 .and.
-     &          transmax(j) <= 0.9999999) then
+              if (transmax(j) .gt. 0.001 .and. 
+     &          transmax(j) .le. 0.9999999) then
                 fluxtopinit(j) = fluxtop(j,ibox)
               tauir(j) = tau(j,ibox) *rec2p13
@@ -829 +829 @@
             do icycle=1,2
               do j=1,npoints  
-                if (tau(j,ibox) > (tauchk            )) then
-                if (transmax(j) > 0.001 .and.
-     &            transmax(j) <= 0.9999999) then
+                if (tau(j,ibox) .gt. (tauchk            )) then 
+                if (transmax(j) .gt. 0.001 .and. 
+     &            transmax(j) .le. 0.9999999) then
                   emcld(j,ibox) = 1. - exp(-1. * tauir(j)  )
                   fluxtop(j,ibox) = fluxtopinit(j) -   
@@ -839 +839 @@
                   tb(j,ibox)= 1307.27
      &              / (log(1. + (1./fluxtop(j,ibox))))
-                  if (tb(j,ibox) > 260.) then
+                  if (tb(j,ibox) .gt. 260.) then
                   tauir(j) = tau(j,ibox) / 2.56
                   end if                   
@@ -850 +850 @@
         
           do j=1,npoints
-            if (tau(j,ibox) > (tauchk            )) then
+            if (tau(j,ibox) .gt. (tauchk            )) then 
                 !cloudy box 
                 !NOTE: tb is the cloud-top temperature not infrared brightness temperature 
                 !at this point in the code
                 tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
-                if (top_height==1.and.tauir(j)<taumin(j)) then
+                if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
                          tb(j,ibox) = attrop(j) - 5. 
                    tau(j,ibox) = 2.13*taumin(j)
@@ -866 +866 @@
         enddo ! ibox
 
-        if (ncolprint/=0) then
+        if (ncolprint.ne.0) then
 
           do j=1,npoints,1000
@@ -925 +925 @@
 
       !compute cloud top pressure
-      do ibox=1,ncol
+      do 30 ibox=1,ncol
         !segregate according to optical thickness
-        if (top_height == 1 .or. top_height == 3) then
+        if (top_height .eq. 1 .or. top_height .eq. 3) then  
           !find level whose temperature
           !most closely matches brightness temperature
@@ -933 +933 @@
             nmatch(j)=0
           enddo
-          do k1=1,nlev-1
-            if (top_height_direction == 2) then
+          do 29 k1=1,nlev-1
+            if (top_height_direction .eq. 2) then
               ilev = nlev - k1 
             else
@@ -941 +941 @@
             !cdir nodep
             do j=1,npoints 
-             if (ilev >= itrop(j)) then
-              if ((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
+             if (ilev .ge. itrop(j)) then
+              if ((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 
                 nmatch(j)=nmatch(j)+1
                 match(j,nmatch(j))=ilev
@@ -951 +951 @@
              end if                         
             enddo
-      END DO
+29        continue
 
           do j=1,npoints 
-            if (nmatch(j) >= 1) then
+            if (nmatch(j) .ge. 1) then
               k1 = match(j,nmatch(j))
               k2 = k1 + 1
@@ -962 +962 @@
               logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
               ptop(j,ibox) = exp(logp)
-              if(abs(pfull(j,k1)-ptop(j,ibox)) <
+              if(abs(pfull(j,k1)-ptop(j,ibox)) .lt.
      &            abs(pfull(j,k2)-ptop(j,ibox))) then
                  levmatch(j,ibox)=k1
@@ -969 +969 @@
               end if   
             else
-              if (tb(j,ibox) <= attrop(j)) then
+              if (tb(j,ibox) .le. attrop(j)) then
                 ptop(j,ibox)=ptrop(j)
                 levmatch(j,ibox)=itrop(j)
               end if
-              if (tb(j,ibox) >= atmax(j)) then
+              if (tb(j,ibox) .ge. atmax(j)) then
                 ptop(j,ibox)=pfull(j,nlev)
                 levmatch(j,ibox)=nlev
@@ -987 +987 @@
           do ilev=1,nlev
             do j=1,npoints     
-              if ((ptop(j,ibox) == 0. )
-     &           .and.(frac_out(j,ibox,ilev) /= 0)) then
+              if ((ptop(j,ibox) .eq. 0. )
+     &           .and.(frac_out(j,ibox,ilev) .ne. 0)) then
                 ptop(j,ibox)=phalf(j,ilev)
               levmatch(j,ibox)=ilev
               end if
-            END DO
-          END DO
+            end do
+          end do
         end if                            
           
         do j=1,npoints
-          if (tau(j,ibox) <= (tauchk            )) then
+          if (tau(j,ibox) .le. (tauchk            )) then
             ptop(j,ibox)=0.
             levmatch(j,ibox)=0      
@@ -1003 +1003 @@
         enddo
 
-      END DO
+30    continue
               
 !
@@ -1032 +1032 @@
 
       !reset frequencies
-      do ilev=1,7
+      do 38 ilev=1,7
       do 38 ilev2=1,7
         do j=1,npoints ! 
-             if (sunlit(j)==1 .or. top_height == 3) then
+             if (sunlit(j).eq.1 .or. top_height .eq. 3) then 
                 fq_isccp(j,ilev,ilev2)= 0.
              else
@@ -1042 +1042 @@
         enddo
 38    continue
-      END DO
 
       !reset variables need for averaging cloud properties
       do j=1,npoints 
-        if (sunlit(j)==1 .or. top_height == 3) then
+        if (sunlit(j).eq.1 .or. top_height .eq. 3) then 
              totalcldarea(j) = 0.
              meanalbedocld(j) = 0.
@@ -1061 +1060 @@
       boxarea = 1./real(ncol)
      
-      do ibox=1,ncol
+      do 39 ibox=1,ncol
         do j=1,npoints 
 
-          if (tau(j,ibox) > (tauchk            )
-     &      .and. ptop(j,ibox) > 0.) then
+          if (tau(j,ibox) .gt. (tauchk            )
+     &      .and. ptop(j,ibox) .gt. 0.) then
               box_cloudy(j,ibox)=.true.
           endif
@@ -1071 +1070 @@
           if (box_cloudy(j,ibox)) then
 
-              if (sunlit(j)==1 .or. top_height == 3) then
+              if (sunlit(j).eq.1 .or. top_height .eq. 3) then
 
                 boxtau(j,ibox) = tau(j,ibox)
 
-                if (tau(j,ibox) >= isccp_taumin) then
+                if (tau(j,ibox) .ge. isccp_taumin) then
                    totalcldarea(j) = totalcldarea(j) + boxarea
                 
@@ -1092 +1091 @@
           endif
 
-          if (sunlit(j)==1 .or. top_height == 3) then
+          if (sunlit(j).eq.1 .or. top_height .eq. 3) then 
 
            if (box_cloudy(j,ibox)) then
@@ -1102 +1101 @@
               boxptop(j,ibox) = ptop(j,ibox)
     
-              if (tau(j,ibox) >= isccp_taumin) then
+              if (tau(j,ibox) .ge. isccp_taumin) then
                 meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
               end if            
@@ -1111 +1110 @@
 
               !determine optical depth category
-              if (tau(j,ibox) < isccp_taumin) then
+              if (tau(j,ibox) .lt. isccp_taumin) then
                   itau(j)=1
-              else if (tau(j,ibox) >= isccp_taumin
+              else if (tau(j,ibox) .ge. isccp_taumin
      &                                    
-     &          .and. tau(j,ibox) < 1.3) then
+     &          .and. tau(j,ibox) .lt. 1.3) then
                 itau(j)=2
-              else if (tau(j,ibox) >= 1.3
-     &          .and. tau(j,ibox) < 3.6) then
+              else if (tau(j,ibox) .ge. 1.3 
+     &          .and. tau(j,ibox) .lt. 3.6) then
                 itau(j)=3
-              else if (tau(j,ibox) >= 3.6
-     &          .and. tau(j,ibox) < 9.4) then
+              else if (tau(j,ibox) .ge. 3.6 
+     &          .and. tau(j,ibox) .lt. 9.4) then
                   itau(j)=4
-              else if (tau(j,ibox) >= 9.4
-     &          .and. tau(j,ibox) < 23.) then
+              else if (tau(j,ibox) .ge. 9.4 
+     &          .and. tau(j,ibox) .lt. 23.) then
                   itau(j)=5
-              else if (tau(j,ibox) >= 23.
-     &          .and. tau(j,ibox) < 60.) then
+              else if (tau(j,ibox) .ge. 23. 
+     &          .and. tau(j,ibox) .lt. 60.) then
                   itau(j)=6
-              else if (tau(j,ibox) >= 60.) then
+              else if (tau(j,ibox) .ge. 60.) then
                   itau(j)=7
               end if
 
               !determine cloud top pressure category
-              if (    ptop(j,ibox) > 0.
-     &          .and.ptop(j,ibox) < 180.) then
+              if (    ptop(j,ibox) .gt. 0.  
+     &          .and.ptop(j,ibox) .lt. 180.) then
                   ipres(j)=1
-              else if(ptop(j,ibox) >= 180.
-     &          .and.ptop(j,ibox) < 310.) then
+              else if(ptop(j,ibox) .ge. 180.
+     &          .and.ptop(j,ibox) .lt. 310.) then
                   ipres(j)=2
-              else if(ptop(j,ibox) >= 310.
-     &          .and.ptop(j,ibox) < 440.) then
+              else if(ptop(j,ibox) .ge. 310.
+     &          .and.ptop(j,ibox) .lt. 440.) then
                   ipres(j)=3
-              else if(ptop(j,ibox) >= 440.
-     &          .and.ptop(j,ibox) < 560.) then
+              else if(ptop(j,ibox) .ge. 440.
+     &          .and.ptop(j,ibox) .lt. 560.) then
                   ipres(j)=4
-              else if(ptop(j,ibox) >= 560.
-     &          .and.ptop(j,ibox) < 680.) then
+              else if(ptop(j,ibox) .ge. 560.
+     &          .and.ptop(j,ibox) .lt. 680.) then
                   ipres(j)=5
-              else if(ptop(j,ibox) >= 680.
-     &          .and.ptop(j,ibox) < 800.) then
+              else if(ptop(j,ibox) .ge. 680.
+     &          .and.ptop(j,ibox) .lt. 800.) then
                   ipres(j)=6
-              else if(ptop(j,ibox) >= 800.) then
+              else if(ptop(j,ibox) .ge. 800.) then
                   ipres(j)=7
               end if 
 
               !update frequencies
-              if(ipres(j) > 0.and.itau(j) > 0) then
+              if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
               fq_isccp(j,itau(j),ipres(j))=
      &          fq_isccp(j,itau(j),ipres(j))+ boxarea
@@ -1167 +1166 @@
                        
         enddo ! j
-      END DO
+39    continue
       
       !compute mean cloud properties
       do j=1,npoints 
-        if (totalcldarea(j) > 0.) then
+        if (totalcldarea(j) .gt. 0.) then
           ! code above guarantees that totalcldarea > 0 
           ! only if sunlit .eq. 1 .or. top_height = 3 
@@ -1194 +1193 @@
 !     OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
 !
-      if (debugcol/=0) then
+      if (debugcol.ne.0) then
 !     
          do j=1,npoints,debugcol
@@ -1208 +1207 @@
               do ibox=1,ncol
                    acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
-                   if (levmatch(j,ibox) == ilev)
+                   if (levmatch(j,ibox) .eq. ilev) 
      &                 acc(ilev,ibox)=acc(ilev,ibox)+1
               enddo
@@ -1228 +1227 @@
      &           (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev) 
      &           ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev) 
-             END DO
+             end do
              close(9)
 
-             if (ncolprint/=0) then
+             if (ncolprint.ne.0) then
                write(6,'(a1)') ' '
                     write(6,'(a2,1X,5(a7,1X),a50)') 

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/lidar_simulator.F90

@@ -247 +247 @@
 !------------------------------------------------------------
 
-      if ( npart /= 4 ) then
+      if ( npart .ne. 4 ) then
         print *,'Error in lidar_simulator, npart should be 4, not',npart
         stop
@@ -267 +267 @@
          polpart(INDX_LSLIQ,5) =  0.0626
 !*     LS Ice coefficients: 
-      if (ice_type==0) then
+      if (ice_type.eq.0) then     
          polpart(INDX_LSICE,1) = -1.0176e-8
          polpart(INDX_LSICE,2) =  1.7615e-6
@@ -275 +275 @@
       endif
 !*     LS Ice NS coefficients: 
-      if (ice_type==1) then
+      if (ice_type.eq.1) then 
          polpart(INDX_LSICE,1) = 1.3615e-8
          polpart(INDX_LSICE,2) = -2.04206e-6
@@ -289 +289 @@
          polpart(INDX_CVLIQ,5) =  0.0626
 !*     CONV Ice coefficients: 
-      if (ice_type==0) then
+      if (ice_type.eq.0) then 
          polpart(INDX_CVICE,1) = -1.0176e-8
          polpart(INDX_CVICE,2) =  1.7615e-6
@@ -296 +296 @@
          polpart(INDX_CVICE,5) =  0.0460
       endif
-      if (ice_type==1) then
+      if (ice_type.eq.1) then
          polpart(INDX_CVICE,1) = 1.3615e-8
          polpart(INDX_CVICE,2) = -2.04206e-6
@@ -342 +342 @@
 ! polynomes kp_lidar derived from Mie theory:
       do i = 1, npart
-       where ( rad_part(:,:,i)>0.0)
+       where ( rad_part(:,:,i).gt.0.0)
          kp_part(:,:,i) = &
             polpart(i,1)*(rad_part(:,:,i)*1e6)**4 &
@@ -362 +362 @@
 ! alpha of particles in each subcolumn:
       do i = 1, npart
-        where ( rad_part(:,:,i)>0.0)
+        where ( rad_part(:,:,i).gt.0.0)
           alpha_part(:,:,i) = 3.0/4.0 * Qscat &
                  * rhoair(:,:) * qpart(:,:,i) &
@@ -378 +378 @@
 !     opt. thick of each layer
       tau_mol(:,1:nlev) = alpha_mol(:,1:nlev) &
-   *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
+         & *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
 !     opt. thick from TOA
       DO k = nlev-1, 1, -1
@@ -390 +390 @@
 !       opt. thick of each layer
         tau_part(:,:,i) = tau_part(:,:,i) &
-   * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
+           & * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
 !       opt. thick from TOA
         DO k = nlev-1, 1, -1 
@@ -400 +400 @@
 !      Upper layer 
        pmol(:,nlev) = beta_mol(:,nlev) / (2.*tau_mol(:,nlev)) &
-   * (1.-exp(-2.0*tau_mol(:,nlev)))
+            & * (1.-exp(-2.0*tau_mol(:,nlev)))
 !      Other layers
        DO k= nlev-1, 1, -1
         tau_mol_lay(:) = tau_mol(:,k)-tau_mol(:,k+1) ! opt. thick. of layer k
-        WHERE (tau_mol_lay(:)>0.)
+        WHERE (tau_mol_lay(:).GT.0.)
           pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1)) / (2.*tau_mol_lay(:)) &
-   * (1.-exp(-2.0*tau_mol_lay(:)))
+            & * (1.-exp(-2.0*tau_mol_lay(:)))
         ELSEWHERE
 !         This must never happend, but just in case, to avoid div. by 0
@@ -429 +429 @@
 !     Upper layer 
       pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
-   * (1.-exp(-2.0*tautot(:,nlev)))
+            & * (1.-exp(-2.0*tautot(:,nlev)))
 
 !     Other layers
       DO k= nlev-1, 1, -1
           tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
-        WHERE (tautot_lay(:)>0.)
+        WHERE (tautot_lay(:).GT.0.)
           pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
-   * (1.-EXP(-2.0*tautot_lay(:)))
+               & * (1.-EXP(-2.0*tautot_lay(:)))
         ELSEWHERE
 !         This must never happend, but just in case, to avoid div. by 0
@@ -468 +468 @@
 !     Upper layer
       pnorm_ice(:,nlev) = betatot_ice(:,nlev) / (2.*tautot_ice(:,nlev)) &
-   * (1.-exp(-2.0*tautot_ice(:,nlev)))
+            & * (1.-exp(-2.0*tautot_ice(:,nlev)))
 
       DO k= nlev-1, 1, -1
           tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1) 
-        WHERE (tautot_lay_ice(:)>0.)
+        WHERE (tautot_lay_ice(:).GT.0.)
          pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))/(2.*tautot_lay_ice(:)) &
-   * (1.-EXP(-2.0*tautot_lay_ice(:)))
+               & * (1.-EXP(-2.0*tautot_lay_ice(:)))
         ELSEWHERE
          pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))
@@ -483 +483 @@
 !     Upper layer
       pnorm_liq(:,nlev) = betatot_liq(:,nlev) / (2.*tautot_liq(:,nlev)) &
-   * (1.-exp(-2.0*tautot_liq(:,nlev)))
+            & * (1.-exp(-2.0*tautot_liq(:,nlev)))
 
       DO k= nlev-1, 1, -1
           tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1) 
-        WHERE (tautot_lay_liq(:)>0.)
+        WHERE (tautot_lay_liq(:).GT.0.)
           pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))/(2.*tautot_lay_liq(:)) &
-   * (1.-EXP(-2.0*tautot_lay_liq(:)))
+               & * (1.-EXP(-2.0*tautot_lay_liq(:)))
         ELSEWHERE
           pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))
@@ -510 +510 @@
 !     Upper layer 
       beta_perp_ice(:,nlev) = pnorm_perp_ice(:,nlev) * (2.*tautot_ice(:,nlev)) &
-   / (1.-exp(-2.0*tautot_ice(:,nlev)))
+            & / (1.-exp(-2.0*tautot_ice(:,nlev)))
 
       DO k= nlev-1, 1, -1
         tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
-        WHERE (tautot_lay_ice(:)>0.)
+        WHERE (tautot_lay_ice(:).GT.0.)
          beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
-   / (1.-exp(-2.0*tautot_lay_ice(:)))
+            & / (1.-exp(-2.0*tautot_lay_ice(:)))
 
         ELSEWHERE
@@ -526 +526 @@
 !     Upper layer 
       beta_perp_liq(:,nlev) = pnorm_perp_liq(:,nlev) * (2.*tautot_liq(:,nlev)) &
-   / (1.-exp(-2.0*tautot_liq(:,nlev)))
+            & / (1.-exp(-2.0*tautot_liq(:,nlev)))
 
       DO k= nlev-1, 1, -1
           tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1) 
-        WHERE (tautot_lay_liq(:)>0.)
+        WHERE (tautot_lay_liq(:).GT.0.)
          beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ max(seuil,EXP(-2.0*tautot_liq(:,k+1))) &
-   * (2.*tautot_lay_liq(:)) / (1.-exp(-2.0*tautot_lay_liq(:)))
+            & * (2.*tautot_lay_liq(:)) / (1.-exp(-2.0*tautot_lay_liq(:)))
 
         ELSEWHERE
@@ -547 +547 @@
 ! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
 !     Upper layer 
-    WHERE(tautot(:,nlev)>0)
+    WHERE(tautot(:,nlev).GT.0)
           pnorm_perp_tot(:,nlev) = &
               (beta_perp_ice(:,nlev)+beta_perp_liq(:,nlev)-(beta_mol(:,nlev)/(1+1/0.0284))) / (2.*tautot(:,nlev)) &
-   * (1.-exp(-2.0*tautot(:,nlev)))
+              & * (1.-exp(-2.0*tautot(:,nlev)))
     ELSEWHERE
     pnorm_perp_tot(:,nlev) = 0.
@@ -563 +563 @@
           ! We remove one contribution using 
           ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following equations:
-            WHERE (pnorm(:,k)==0)
+            WHERE (pnorm(:,k).eq.0)
                   pnorm_perp_tot(:,k)=0.
                   ELSEWHERE
-                    WHERE (tautot_lay(:)>0.)
+                    WHERE (tautot_lay(:).GT.0.)
                       pnorm_perp_tot(:,k) = &
                           (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
                           EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
-   * (1.-EXP(-2.0*tautot_lay(:)))
+                          & * (1.-EXP(-2.0*tautot_lay(:)))
                     ELSEWHERE
           !         This must never happen, but just in case, to avoid div. by 0
@@ -690 +690 @@
 ! Lum_norm=f(tetaS,tau_cloud) derived from adding-doubling calculations
 !        valid ONLY ABOVE OCEAN (albedo_sfce=5%)
-!        valid only in one viewing direction (theta_v=30°, phi_s-phi_v=320°)
+!        valid only in one viewing direction (theta_v=30�, phi_s-phi_v=320�)
 !        based on adding-doubling radiative transfer computation
 !        for tau values (0 to 100) and for tetas values (0 to 80)
 !        for 2 scattering phase functions: liquid spherical, ice non spherical
 
-    IF ( nrefl> ntetas ) THEN
+    IF ( nrefl.GT. ntetas ) THEN
         PRINT *,'Error in lidar_simulator, nrefl should be less then ',ntetas,' not',nrefl
         STOP
@@ -711 +711 @@
 !
 ! relative fraction of the opt. thick due to liquid or ice clouds
-    WHERE (tautot_S(:) > 0.)
+    WHERE (tautot_S(:) .GT. 0.)
         frac_taucol_liq(:) = tautot_S_liq(:) / tautot_S(:)
         frac_taucol_ice(:) = tautot_S_ice(:) / tautot_S(:)
@@ -733 +733 @@
     DO it=1,ntetas
       DO ny=1,nbtau-1
-        WHERE (tautot_S(:)>=tau(ny).AND.tautot_S(:)<=tau(ny+1))
+        WHERE (tautot_S(:).GE.tau(ny).AND.tautot_S(:).LE.tau(ny+1))
             rlumA_mod(:,it) = aA(it,ny)*tautot_S(:) + bA(it,ny)
             rlumB_mod(:,it) = aB(it,ny)*tautot_S(:) + bB(it,ny)

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90

@@ -33 +33 @@
          IRNGT (1) = 1
          Return
+      Case Default
+         Continue
       End Select
 !
@@ -231 +233 @@
          IRNGT (1) = 1
          Return
+      Case Default
+         Continue
       End Select
 !

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/math_lib.F90

@@ -34 +34 @@
 
 ! ----- INPUTS -----
-  real(kind=8), intent(in) :: x
+  real*8, intent(in) :: x
   
 ! ----- OUTPUTS -----
-  real(kind=8) :: gamma
+  real*8 :: gamma
 
 ! ----- INTERNAL -----  
-  real(kind=8) :: pi,ga,z,r,gr
-  real(kind=8) :: g(26)
+  real*8 :: pi,ga,z,r,gr
+  real*8 :: g(26)
   integer :: k,m1,m
        
@@ -124 +124 @@
 
 ! ----- INPUTS -----  
-  real(kind=8), intent(in), dimension(:) :: f,s
+  real*8, intent(in), dimension(:) :: f,s  
   integer, intent(in) :: i1, i2
 
 ! ---- OUTPUTS -----
-  real(kind=8) :: path_integral
+  real*8 :: path_integral  
   
 ! ----- INTERNAL -----    
-  real(kind=8) :: sumo, deltah, val
-  integer(kind=4) :: nelm, j
-  integer(kind=4), dimension(i2-i1+1) :: idx
-  real(kind=8), dimension(i2-i1+1) :: f_rev, s_rev
+  real*8 :: sumo, deltah, val
+  integer*4 :: nelm, j
+  integer*4, dimension(i2-i1+1) :: idx
+  real*8, dimension(i2-i1+1) :: f_rev, s_rev
 
   nelm = i2-i1+1
@@ -273 +273 @@
        exit
     end if
-  END DO
+  end do
   
   if (lerror) then
@@ -316 +316 @@
     end if
     ilo = ilo + 1
-  END DO
+  end do
 
   ilo = max ( 2, ilo )
@@ -326 +326 @@
     end if
     ihi = ihi - 1
-  END DO
+  end do
   
   ihi = min ( ihi, ntab - 1 )
@@ -374 +374 @@
     syl = x2
  
-  END DO
+  end do
  
   result = sum1 &

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp.F90

@@ -179 +179 @@
    minp = minval(gbx%psfc)
    maxp = maxval(gbx%psfc)
-   if (Npoints > 1) seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
+   if (Npoints .gt. 1) seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
    ! Below it's how it was done in the original implementation of the ISCCP simulator. 
    ! The one above is better for offline data, when you may have packed data 
@@ -414 +414 @@
                 if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
                 if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
-                if (sgx%prec_frac(j,i,Nlevels+1-k) == 1) prec_ls(j,k)=prec_ls(j,k)+1.
-                if (sgx%prec_frac(j,i,Nlevels+1-k) == 2) prec_cv(j,k)=prec_cv(j,k)+1.
-                if (sgx%prec_frac(j,i,Nlevels+1-k) == 3) then
+                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
+                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
+                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 3) then
                     prec_cv(j,k)=prec_cv(j,k)+1.
                     prec_ls(j,k)=prec_ls(j,k)+1.
@@ -501 +501 @@
             do j=1,Npoints
                 !--------- Clouds -------
-                if (frac_ls(j,k) /= 0.) then
+                if (frac_ls(j,k) .ne. 0.) then
                     sghydro%mr_hydro(j,:,k,I_LSCLIQ) = sghydro%mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
                     sghydro%mr_hydro(j,:,k,I_LSCICE) = sghydro%mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
                 endif
-                if (frac_cv(j,k) /= 0.) then
+                if (frac_cv(j,k) .ne. 0.) then
                     sghydro%mr_hydro(j,:,k,I_CVCLIQ) = sghydro%mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
                     sghydro%mr_hydro(j,:,k,I_CVCICE) = sghydro%mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
@@ -511 +511 @@
                 !--------- Precip -------
                 if (gbx%use_precipitation_fluxes) then
-                    if (prec_ls(j,k) /= 0.) then
+                    if (prec_ls(j,k) .ne. 0.) then
                         gbx%rain_ls(j,k) = gbx%rain_ls(j,k)/prec_ls(j,k)
                         gbx%snow_ls(j,k) = gbx%snow_ls(j,k)/prec_ls(j,k)
                         gbx%grpl_ls(j,k) = gbx%grpl_ls(j,k)/prec_ls(j,k)
                     endif
-                    if (prec_cv(j,k) /= 0.) then
+                    if (prec_cv(j,k) .ne. 0.) then
                         gbx%rain_cv(j,k) = gbx%rain_cv(j,k)/prec_cv(j,k)
                         gbx%snow_cv(j,k) = gbx%snow_cv(j,k)/prec_cv(j,k)
                     endif
                 else
-                    if (prec_ls(j,k) /= 0.) then
+                    if (prec_ls(j,k) .ne. 0.) then
                         sghydro%mr_hydro(j,:,k,I_LSRAIN) = sghydro%mr_hydro(j,:,k,I_LSRAIN)/prec_ls(j,k)
                         sghydro%mr_hydro(j,:,k,I_LSSNOW) = sghydro%mr_hydro(j,:,k,I_LSSNOW)/prec_ls(j,k)
                         sghydro%mr_hydro(j,:,k,I_LSGRPL) = sghydro%mr_hydro(j,:,k,I_LSGRPL)/prec_ls(j,k)
                     endif
-                    if (prec_cv(j,k) /= 0.) then
+                    if (prec_cv(j,k) .ne. 0.) then
                         sghydro%mr_hydro(j,:,k,I_CVRAIN) = sghydro%mr_hydro(j,:,k,I_CVRAIN)/prec_cv(j,k)
                         sghydro%mr_hydro(j,:,k,I_CVSNOW) = sghydro%mr_hydro(j,:,k,I_CVSNOW)/prec_cv(j,k)

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_modis_simulator.F90

@@ -176 +176 @@
               opticalThickness(i, j, k) = 0.   
             end if 
-          END DO
-        END DO
-      END DO
+          end do 
+        end do
+      end do
 
       !
@@ -197 +197 @@
           do i = 1, nSunlit
             if(subCols%frac_out(sunlit(i), j, k) == I_CVC) opticalThickness(i, j, k) = gridBox%dtau_c(sunlit(i), k)
-          END DO
-        END DO
-      END DO
+          end do 
+        end do
+      end do
 
       !
@@ -220 +220 @@
                                 retrievedPhase(i, :), retrievedCloudTopPressure(i, :),      & 
                                 retrievedTau(i, :), retrievedSize(i, :))
-     END DO
+     end do
      
       ! DJS2015: Call L3 modis simulator used by cospv2.0

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_radar.F90

@@ -53 +53 @@
 
         real undef
-        real(kind=8), dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
+        real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
             t_matrix,rh_matrix
-        real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
-        real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
-        real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(inout) :: Np_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
+        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: Np_matrix
 
         ! ----- OUTPUTS -----
-        real(kind=8), dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
+        real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
             g_to_vol,dBZe,a_to_vol
         ! ----- OPTIONAL -----
-        real(kind=8), optional, dimension(nprof,ngate) :: &
+        real*8, optional, dimension(nprof,ngate) :: &
             g_to_vol_in,g_to_vol_out
      end subroutine radar_simulator
@@ -86 +86 @@
   nsizes            ! num of discrete drop sizes
 
-  real(kind=8), dimension(:,:), allocatable :: &
+  real*8, dimension(:,:), allocatable :: &
   g_to_vol ! integrated atten due to gases, r>v (dB)
 
-  real(kind=8), dimension(:,:), allocatable :: &
+  real*8, dimension(:,:), allocatable :: &
   Ze_non, &         ! radar reflectivity withOUT attenuation (dBZ)
   Ze_ray, &         ! Rayleigh reflectivity (dBZ)
@@ -100 +100 @@
   rh_matrix                     !relative humidity (%)
 
-  real(kind=8), dimension(:,:,:), allocatable :: &
+  real*8, dimension(:,:,:), allocatable :: &
   hm_matrix, &          ! hydrometeor mixing ratio (g kg^-1)
   re_matrix, &          ! effective radius (microns).   Optional. 0 ==> use Np_matrix or defaults

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_types.F90

@@ -1110 +1110 @@
 
 
-    if(y%Nhydro/=N_HYDRO) then
+    if(y%Nhydro.ne.N_HYDRO) then
 
         write(*,*) 'Number of hydrometeor input to subroutine', &

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_llnl_stats.F90

@@ -119 +119 @@
        do j=Nlevels,1,-1 !top->surf
         sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
-        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
+        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
             lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
             flag_cld=1
@@ -129 +129 @@
         endif
        enddo !levels
-       if (flag_cld == 1) tcc(pr)=tcc(pr)+1.
+       if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1.
      enddo !columns
    enddo !points

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_lmd_ipsl_stats.F90

@@ -148 +148 @@
       do ic = 1, ncol
         pnorm_c = pnorm(:,ic,:)
-        where ((pnorm_c<xmax) .and. (pmol<xmax) .and. (pmol> 0.0 ))
+        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
             x3d_c = pnorm_c/pmol
         elsewhere
@@ -247 +247 @@
 ! c 0- Initializations
 ! c -------------------------------------------------------
-      if ( Nbins < 6) return
+      if ( Nbins .lt. 6) return
 
       srbval(1) =  S_att
@@ -275 +275 @@
                do i = 1, Npoints
                   if (x(i,k,j) /= undef) then
-                     if ((x(i,k,j)>srbval_ext(ib-1)).and.(x(i,k,j)<=srbval_ext(ib))) &
+                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
                           cfad(i,ib,j) = cfad(i,ib,j) + 1.0
                   else 
@@ -285 +285 @@
       enddo
 
-      where (cfad /= undef)  cfad = cfad / float(Ncolumns)
+      where (cfad .ne. undef)  cfad = cfad / float(Ncolumns)
 
 ! c -------------------------------------------------------
@@ -373 +373 @@
 ! ---------------------------------------------------------------
 
-      if ( Ncat /= 4 ) then
+      if ( Ncat .ne. 4 ) then
          print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat
          stop
@@ -423 +423 @@
 
 ! cloud detection at subgrid-scale:
-         where ( (x(:,:,k)>S_cld) .and. (x(:,:,k)/= undef) )
+         where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) )
            cldy(:,:,k)=1.0
          elsewhere
@@ -430 +430 @@
 
 ! number of usefull sub-columns:
-         where ( (x(:,:,k)>S_att) .and. (x(:,:,k)/= undef)  )
+         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
            srok(:,:,k)=1.0
          elsewhere
@@ -513 +513 @@
           ! 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)>0.)then
+          if(srok(ip,ic,k).gt.0.)then
           do itemp=1,Ntemp
-            if( (tmp(ip,k)>=tempmod(itemp)).and.(tmp(ip,k)<tempmod(itemp+1)) )then
+            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
               lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
             endif
@@ -521 +521 @@
           endif
 
-          if(cldy(ip,ic,k)==1.)then
+          if(cldy(ip,ic,k).eq.1.)then
           do itemp=1,Ntemp
-            if( (tmp(ip,k)>=tempmod(itemp)).and.(tmp(ip,k)<tempmod(itemp+1)) )then
+            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
               lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
             endif
@@ -532 +532 @@
           iz=1
           p1 = pplay(ip,k)
-          if ( p1>0. .and. p1<(440.*100.)) then ! high clouds
+          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
             iz=3
-          else if(p1>=(440.*100.) .and. p1<(680.*100.)) then  ! mid clouds
+          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
             iz=2
          endif
@@ -554 +554 @@
 ! -- grid-box 3D cloud fraction
 
-      where ( nsub(:,:)>0.0 )
+      where ( nsub(:,:).gt.0.0 )
          lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
       elsewhere
@@ -573 +573 @@
        enddo
       enddo
-      where ( nsublayer(:,:)>0.0 )
+      where ( nsublayer(:,:).gt.0.0 )
          cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
       elsewhere
@@ -593 +593 @@
 
 ! Avoid zero values
-        if( (cldy(i,ncol,nlev)==1.) .and. (ATBperp(i,ncol,nlev)>0.) )then
+        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
 ! Computation of the ATBperp along the phase discrimination line
            ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
@@ -604 +604 @@
 !____________________________________________________________________________________________________
 !
-           if( (ATBperp(i,ncol,nlev)-ATBperp_tmp)>=0. )then   ! Ice clouds
+           if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
              ! ICE with temperature above 273,15°K = Liquid (false ice)
-            if(tmp(i,nlev)>273.15)then                ! Temperature above 273,15 K
+            if(tmp(i,nlev).gt.273.15)then                ! Temperature above 273,15 K
               ! Liquid: False ice corrected by the temperature to Liquid
                lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.   ! false ice detection ==> added to Liquid
@@ -613 +613 @@
                                                     ! to classify the phase cloud
                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,2) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,2) = 1.
                 else                                                    ! low cloud
@@ -621 +621 @@
                 endif
                    cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,5) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,5) = 1.
                 else                                                    ! low cloud
@@ -634 +634 @@
               tmpi(i,ncol,nlev)=tmp(i,nlev)
                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,1) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,1) = 1.
                 else                                                    ! low cloud
@@ -651 +651 @@
              else                                        ! Liquid clouds
               ! Liquid with temperature above 231,15°K
-            if(tmp(i,nlev)>231.15)then
+            if(tmp(i,nlev).gt.231.15)then 
                lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
                tmpl(i,ncol,nlev)=tmp(i,nlev)
                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,2) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,2) = 1.
                 else                                                    ! low cloud
@@ -670 +670 @@
                                                     ! to classify the phase cloud
                    cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,4) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,4) = 1.
                 else                                                    ! low cloud
@@ -678 +678 @@
                 endif
                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,1) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,1) = 1.
                 else                                                    ! low cloud
@@ -702 +702 @@
          p1 = pplay(i,nlev)
 
-        if( (cldy(i,ncol,nlev)==1.) .and. (ATBperp(i,ncol,nlev)>0.) )then
+        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
 ! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 
 !                                  + ATB*epsilon50 + zeta50
@@ -715 +715 @@
 !
             ! ICE with temperature above 273,15°K = Liquid (false ice)
-          if( (ATBperp(i,ncol,nlev)-ATBperp_tmp)>=0. )then   ! Ice clouds
-            if(tmp(i,nlev)>273.15)then
+          if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
+            if(tmp(i,nlev).gt.273.15)then 
                lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.  ! false ice ==> liq
                tmpl(i,ncol,nlev)=tmp(i,nlev)
@@ -722 +722 @@
 
                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
-               if ( p1>0. .and. p1<(440.*100.)) then              ! high cloud
+               if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
                    cldlayphase(i,ncol,3,2) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,2) = 1.
                 else                                                    ! low cloud
@@ -731 +731 @@
 
                    cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,5) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,5) = 1.
                 else                                                    ! low cloud
@@ -745 +745 @@
 
                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,1) = 1.
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,1) = 1.
                 else                                                    ! low cloud
@@ -762 +762 @@
           else  
              ! Liquid with temperature above 231,15°K
-            if(tmp(i,nlev)>231.15)then
+            if(tmp(i,nlev).gt.231.15)then 
                lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
                tmpl(i,ncol,nlev)=tmp(i,nlev)
 
                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,2) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,2) = 1.
                 else                                                    ! low cloud
@@ -782 +782 @@
 
                    cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,4) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,4) = 1.
                 else                                                    ! low cloud
@@ -791 +791 @@
 
                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
-                if ( p1>0. .and. p1<(440.*100.)) then             ! high cloud
+                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
                    cldlayphase(i,ncol,3,1) = 1.  
-                else if(p1>=(440.*100.) .and. p1<(680.*100.)) then ! mid cloud
+                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
                    cldlayphase(i,ncol,2,1) = 1.
                 else                                                    ! low cloud
@@ -805 +805 @@
 
            ! Find the level of the highest cloud with SR>30
-            if(x(i,ncol,nlev)>S_cld_att)then     ! SR > 30.
+            if(x(i,ncol,nlev).gt.S_cld_att)then  ! SR > 30.
                 toplvlsat=nlev-1
                 goto 99 
@@ -821 +821 @@
 !____________________________________________________________________________________________________
 !
-if(toplvlsat/=0)then
+if(toplvlsat.ne.0)then          
       do nlev=toplvlsat,1,-1
          p1 = pplay(i,nlev)
-        if(cldy(i,ncol,nlev)==1.)then
+        if(cldy(i,ncol,nlev).eq.1.)then
            lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
            tmpu(i,ncol,nlev)=tmp(i,nlev)
 
                    cldlayphase(i,ncol,4,3) = 1.                         ! tot cloud
-          if ( p1>0. .and. p1<(440.*100.)) then              ! high cloud
+          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
              cldlayphase(i,ncol,3,3) = 1.
-          else if(p1>=(440.*100.) .and. p1<(680.*100.)) then  ! mid cloud
+          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid cloud
              cldlayphase(i,ncol,2,3) = 1.
           else                                                     ! low cloud
@@ -857 +857 @@
 ! of the occurrences
 lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
-WHERE (lidarcldphasetmp(:,:)> 0.)
+WHERE (lidarcldphasetmp(:,:).gt. 0.)
    lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
 ELSEWHERE
@@ -864 +864 @@
 
 ! Compute Phase 3D Cloud Fraction
-     WHERE ( nsub(:,:)>0.0 )
+     WHERE ( nsub(:,:).gt.0.0 )
        lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
        lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
@@ -899 +899 @@
 ! Compute the Ice percentage in cloud = ice/(ice+liq)
 cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
-    WHERE (cldlayerphasetmp(:,:)> 0.)
+    WHERE (cldlayerphasetmp(:,:).gt. 0.)
        cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
     ELSEWHERE
@@ -906 +906 @@
 
     do i=1,Nphase-1
-      WHERE ( cldlayerphasesum(:,:)>0.0 )
+      WHERE ( cldlayerphasesum(:,:).gt.0.0 )
          cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:) 
       ENDWHERE
@@ -917 +917 @@
           checkcldlayerphase2=0.
 
-          if (cldlayerphasesum(i,iz)>0.0 )then
+          if (cldlayerphasesum(i,iz).gt.0.0 )then
              do ic=1,Nphase-3
                 checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic)  
              enddo
              checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase
-             if( (checkcldlayerphase2>0.01).or.(checkcldlayerphase2<-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
+             if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
 
           endif
@@ -930 +930 @@
 
     do i=1,Nphase-1
-      WHERE ( nsublayer(:,:)==0.0 )
+      WHERE ( nsublayer(:,:).eq.0.0 )
          cldlayerphase(:,:,i) = undef
       ENDWHERE
@@ -942 +942 @@
 do i=1,Npoints
 do itemp=1,Ntemp
-if(tmpi(i,ncol,nlev)>0.)then
-      if( (tmpi(i,ncol,nlev)>=tempmod(itemp)).and.(tmpi(i,ncol,nlev)<tempmod(itemp+1)) )then
+if(tmpi(i,ncol,nlev).gt.0.)then
+      if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
         lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1.
       endif
-elseif(tmpl(i,ncol,nlev)>0.)then
-      if( (tmpl(i,ncol,nlev)>=tempmod(itemp)).and.(tmpl(i,ncol,nlev)<tempmod(itemp+1)) )then
+elseif(tmpl(i,ncol,nlev).gt.0.)then
+      if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then
         lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1.
       endif
-elseif(tmpu(i,ncol,nlev)>0.)then
-      if( (tmpu(i,ncol,nlev)>=tempmod(itemp)).and.(tmpu(i,ncol,nlev)<tempmod(itemp+1)) )then
+elseif(tmpu(i,ncol,nlev).gt.0.)then
+      if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then
         lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1.
       endif
@@ -965 +965 @@
 checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
 
-        if(checktemp/=lidarcldtemp(i,itemp,1))then
+        if(checktemp.NE.lidarcldtemp(i,itemp,1))then
           print *, i,itemp
           print *, lidarcldtemp(i,itemp,1:4)
@@ -984 +984 @@
 
 do i=1,4
-  WHERE(lidarcldtempind(:,:)>0.)
+  WHERE(lidarcldtempind(:,:).gt.0.)
      lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
   ELSEWHERE
@@ -1046 +1046 @@
     do k=1,Nlevels
        ! Cloud detection at subgrid-scale:
-       where ( (x(:,:,k) > S_cld) .and. (x(:,:,k) /= undef) )
+       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
           cldy(:,:,k)=1.0
        elsewhere
@@ -1052 +1052 @@
        endwhere
        ! Fully attenuated layer detection at subgrid-scale:
-       where ( (x(:,:,k) > 0.0) .and. (x(:,:,k) < S_att_opaq) .and. (x(:,:,k) /= undef) )
+       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ne. undef) )
           cldyopaq(:,:,k)=1.0
        elsewhere
@@ -1059 +1059 @@
 
        ! Number of useful sub-column layers:
-       where ( (x(:,:,k) > S_att) .and. (x(:,:,k) /= undef) )
+       where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
           srok(:,:,k)=1.0
        elsewhere
@@ -1065 +1065 @@
        endwhere
        ! Number of useful sub-columns layers for z_opaque 3D fraction:
-       where ( (x(:,:,k) > 0.0) .and. (x(:,:,k) /= undef) )
+       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .ne. undef) )
           srokopaq(:,:,k)=1.0
        elsewhere
@@ -1098 +1098 @@
 
      ! Declaring non-opaque cloudy profiles as thin cloud profiles
-           if ( (cldlay(ip,ic,4) == 1.0) .and. (cldlay(ip,ic,1) == 0.0) ) then
+           if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
               cldlay(ip,ic,2)  =  1.0
            endif
@@ -1105 +1105 @@
 
      ! Opaque cloud profiles
-           if ( cldlay(ip,ic,1) == 1.0 ) then
+           if ( cldlay(ip,ic,1) .eq. 1.0 ) then
               zopac = 0.0
               do k=2,Nlevels
      ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables
-                 if ( (cldy(ip,ic,k) == 1.0) .and. (zopac == 0.0) ) then
+                 if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
                     lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
                     cldlay(ip,ic,3)        = vgrid_z(k-1) !z_opaque altitude
@@ -1115 +1115 @@
                     zopac = 1.0
                  endif
-                 if ( cldy(ip,ic,k) == 1.0 ) then
+                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
                     lidarcldtype(ip,k,1)   = lidarcldtype(ip,k,1) + 1.0
                  endif
@@ -1122 +1122 @@
 
      ! Thin cloud profiles
-           if ( cldlay(ip,ic,2) == 1.0 ) then
+           if ( cldlay(ip,ic,2) .eq. 1.0 ) then
               do k=1,Nlevels
      ! Declaring thin cloud fraction for 3D variable
-                 if ( cldy(ip,ic,k) == 1.0 ) then
+                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
                     lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0
                  endif
@@ -1135 +1135 @@
 
     ! 3D cloud types fraction (opaque=1 and thin=2)
-    where ( nsub(:,:) > 0.0 )
+    where ( nsub(:,:) .gt. 0.0 )
        lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
        lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
@@ -1143 +1143 @@
     endwhere
     ! 3D z_opaque fraction (=3)
-    where ( nsubopaq(:,:) > 0.0 )
+    where ( nsubopaq(:,:) .gt. 0.0 )
        lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
     elsewhere
@@ -1152 +1152 @@
     do ip = 1, Npoints
         do k = Nlevels-1, 1, -1
-           if ( lidarcldtype(ip,k,3) /= undef ) then
+           if ( lidarcldtype(ip,k,3) .ne. undef ) then
               lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
            endif
         enddo
     enddo
-    where ( nsubopaq(:,:) == 0.0 )
+    where ( nsubopaq(:,:) .eq. 0.0 )
        lidarcldtype(:,:,4) = undef
     endwhere
@@ -1169 +1169 @@
        enddo
     enddo
-    where (nsublayer(:,:) > 0.0)
+    where (nsublayer(:,:) .gt. 0.0)
        cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
     elsewhere

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_modis_sim.F90

@@ -331 +331 @@
         retrievedTau(i) = R_UNDEF 
       end if
-    END DO
+    end do
     where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill
 
@@ -666 +666 @@
     do ij=2,nbin1+1
        do ik=2,nbin2+1
-          jointHist(ij-1,ik-1)=count(var1 >= bin1(ij-1) .and. var1 < bin1(ij) .and. &
-               var2 >= bin2(ik-1) .and. var2 < bin2(ik))
+          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))        
        enddo
     enddo
@@ -802 +802 @@
         tauMask(:, :, i) = .false.
       end where
-    END DO
+    end do 
 
     do i = 1, numPressureHistogramBins 
@@ -811 +811 @@
         pressureMask(:, :, i) = .false.
       end where
-    END DO
+    end do 
     
     do i = 1, numPressureHistogramBins
@@ -817 +817 @@
         Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = & 
           real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
-      END DO
-    END DO
+      end do 
+    end do 
     
   end subroutine modis_L3_simulator
@@ -851 +851 @@
       end if 
       if(totalTau >= tauLimit) exit
-    END DO
+    end do 
     cloud_top_pressure = totalProduct/totalTau
   end function cloud_top_pressure
@@ -877 +877 @@
       end if 
       if(totalTau >= tauLimit) exit
-    END DO
+    end do 
     weight_by_extinction = totalProduct/totalTau
   end function weight_by_extinction
@@ -1114 +1114 @@
     do i = 1, size(cloudIndicies)
       call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
-    END DO
+    end do 
                     
     call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)  
@@ -1292 +1292 @@
           Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i))
           Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1 - Refl_cumulative(i-1) * Refl(i))
-      END DO
+      end do
       
       Refl_tot = Refl_cumulative(size(Refl))

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90

@@ -38 +38 @@
   
 ! ----- INPUTS -----
-  real(kind=8), intent(in) :: freq,tk
+  real*8, intent(in) :: freq,tk
   
 ! ----- OUTPUTS -----
-  real(kind=8), intent(out) :: n_r, n_i
+  real*8, intent(out) :: n_r, n_i
 
 ! ----- INTERNAL -----    
-  real(kind=8) ld,es,ei,a,ls,sg,tm1,cos1,sin1
-  real(kind=8) e_r,e_i
-  real(kind=8) pi
-  real(kind=8) tc
-  complex(kind=8) e_comp, sq
+  real*8 ld,es,ei,a,ls,sg,tm1,cos1,sin1
+  real*8 e_r,e_i
+  real*8 pi
+  real*8 tc
+  complex*16 e_comp, sq
 
   tc = tk - 273.15
@@ -102 +102 @@
 
 ! ----- INPUTS -----
-  real(kind=8), intent(in) :: freq, t
+  real*8, intent(in) :: freq, t
   
 ! ----- OUTPUTS -----  
-  real(kind=8), intent(out) :: n_r,n_i
+  real*8, intent(out) :: n_r,n_i
 
 ! Parameters:
-  integer(kind=2) :: i,lt1,lt2,nwl,nwlt
+  integer*2 :: i,lt1,lt2,nwl,nwlt
   parameter(nwl=468,nwlt=62)
 
-  real(kind=8) :: alam,cutice,pi,t1,t2,wlmax,wlmin, &
+  real*8 :: alam,cutice,pi,t1,t2,wlmax,wlmin, &
             x,x1,x2,y,y1,y2,ylo,yhi,tk
 
-  real(kind=8) :: &
+  real*8 :: &
        tabim(nwl),tabimt(nwlt,4),tabre(nwl),tabret(nwlt,4),temref(4), &
        wl(nwl),wlt(nwlt)
@@ -519 +519 @@
 
 !   // region from 0.045 microns to 167.0 microns - no temperature depend
+    do i=2,nwl
+      if(alam < wl(i)) continue
+    enddo
     x1=log(wl(i-1))
     x2=log(wl(i))
@@ -536 +539 @@
     if(tk > temref(1)) tk=temref(1)
     if(tk < temref(4)) tk=temref(4)
-    do i=2,4
-      if(tk>=temref(i)) go to 12
-    END DO
+    do 11 i=2,4
+      if(tk.ge.temref(i)) go to 12
+    11 continue
     12 lt1=i
     lt2=i-1
-    do i=2,nwlt
-      if(alam<=wlt(i)) go to 14
-    END DO
+    do 13 i=2,nwlt
+      if(alam.le.wlt(i)) go to 14
+    13 continue
     14 x1=log(wlt(i-1))
     x2=log(wlt(i))
@@ -583 +586 @@
       Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error)
 
-      Integer (kind=2)  Imaxx
+      Integer * 2  Imaxx
       Parameter (Imaxx = 12000)
-      Real (kind=4)     RIMax          ! largest real part of refractive index
+      Real * 4     RIMax          ! largest real part of refractive index
       Parameter (RIMax = 2.5)
-      Real (kind=4)     IRIMax         ! largest imaginary part of refractive index
+      Real * 4     IRIMax         ! largest imaginary part of refractive index
       Parameter (IRIMax = -2)
-      Integer (kind=2)  Itermax
+      Integer * 2  Itermax
       Parameter (Itermax = 12000 * 2.5)
                                 ! must be large enough to cope with the
                                 ! largest possible nmx = x * abs(scm) + 15
                                 ! or nmx =  Dx + 4.05*Dx**(1./3.) + 2.0
-      Integer (kind=2)  Imaxnp
+      Integer * 2  Imaxnp
       Parameter (Imaxnp = 10000)  ! Change this as required
 !     INPUT
-      Real (kind=8)     Dx
-      Complex (kind=8)  SCm
-      Integer (kind=4)  Inp
-      Real (kind=8)     Dqv(Inp)
+      Real * 8     Dx
+      Complex * 16  SCm
+      Integer * 4  Inp
+      Real * 8     Dqv(Inp)
 !     OUTPUT
-      Complex (kind=8)  Xs1(InP)
-      Complex (kind=8)  Xs2(InP)
-      Real (kind=8)     Dqxt
-      Real (kind=8)     Dqsc
-      Real (kind=8)     Dg
-      Real (kind=8)     Dbsc
-      Real (kind=8)     DPh(InP)
-      Integer (kind=4)  Error
+      Complex * 16  Xs1(InP)
+      Complex * 16  Xs2(InP)
+      Real * 8     Dqxt
+      Real * 8     Dqsc
+      Real * 8     Dg
+      Real * 8     Dbsc
+      Real * 8     DPh(InP)
+      Integer * 4  Error
 !     LOCAL
-      Integer (kind=2)  I
-      Integer (kind=2)  NStop
-      Integer (kind=2)  NmX
-      Integer (kind=4)  N    ! N*N > 32767 ie N > 181
-      Integer (kind=4)  Inp2
-      Real (kind=8)     Chi,Chi0,Chi1
-      Real (kind=8)     APsi,APsi0,APsi1
-      Real (kind=8)     Pi0(Imaxnp)
-      Real (kind=8)     Pi1(Imaxnp)
-      Real (kind=8)     Taun(Imaxnp)
-      Real (kind=8)     Psi,Psi0,Psi1
-      Complex (kind=4)  Ir
-      Complex (kind=8) Cm
-      Complex (kind=8) A,ANM1,APB
-      Complex (kind=8) B,BNM1,AMB
-      Complex (kind=8) D(Itermax)
-      Complex (kind=8) Sp(Imaxnp)
-      Complex (kind=8) Sm(Imaxnp)
-      Complex (kind=8) Xi,Xi0,Xi1
-      Complex (kind=8) Y
+      Integer * 2  I
+      Integer * 2  NStop
+      Integer * 2  NmX
+      Integer * 4  N    ! N*N > 32767 ie N > 181
+      Integer * 4  Inp2
+      Real * 8     Chi,Chi0,Chi1
+      Real * 8     APsi,APsi0,APsi1
+      Real * 8     Pi0(Imaxnp)
+      Real * 8     Pi1(Imaxnp)
+      Real * 8     Taun(Imaxnp)
+      Real * 8     Psi,Psi0,Psi1
+      Complex * 8  Ir
+      Complex * 16 Cm
+      Complex * 16 A,ANM1,APB
+      Complex * 16 B,BNM1,AMB
+      Complex * 16 D(Itermax)
+      Complex * 16 Sp(Imaxnp)
+      Complex * 16 Sm(Imaxnp)
+      Complex * 16 Xi,Xi0,Xi1
+      Complex * 16 Y
 !     ACCELERATOR VARIABLES
-      Integer (kind=2)  Tnp1
-      Integer (kind=2)  Tnm1
-      Real (kind=8)     Dn
-      Real (kind=8)     Rnx
-      Real (kind=8)     S(Imaxnp)
-      Real (kind=8)     T(Imaxnp)
-      Real (kind=8)     Turbo
-      Real (kind=8)     A2
-      Complex (kind=8) A1
+      Integer * 2  Tnp1
+      Integer * 2  Tnm1
+      Real * 8     Dn
+      Real * 8     Rnx
+      Real * 8     S(Imaxnp)
+      Real * 8     T(Imaxnp)
+      Real * 8     Turbo
+      Real * 8     A2
+      Complex * 16 A1
       
-      If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then
+      If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
         Error = 1
         Return
@@ -649 +652 @@
       Ir = 1 / Cm
       Y =  Dx * Cm
-      If (Dx<0.02) Then
+      If (Dx.Lt.0.02) Then
          NStop = 2
       Else
-         If (Dx<=8.0) Then
+         If (Dx.Le.8.0) Then
             NStop = Dx + 4.00*Dx**(1./3.) + 2.0
          Else
-            If (Dx< 4200.0) Then
+            If (Dx.Lt. 4200.0) Then
                NStop = Dx + 4.05*Dx**(1./3.) + 2.0
             Else
@@ -663 +666 @@
       End If
       NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
-      If (Nmx > Itermax) then
+      If (Nmx .gt. Itermax) then
           Error = 1
           Return
@@ -706 +709 @@
          Dqxt = Tnp1 *      Dble(A + B)          + Dqxt
          Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
-         If (N>1) then
+         If (N.Gt.1) then
          Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
          End If
@@ -714 +717 @@
          AMB = A2 * (A - B)
          Do I = 1,Inp2
-            If (I>Inp) Then
+            If (I.GT.Inp) Then
                S(I) = -Pi1(I)
             Else
@@ -733 +736 @@
          Xi1 = Dcmplx(APsi1,Chi1)
       End Do
-      If (Dg >0) Dg = 2 * Dg / Dqsc
+      If (Dg .GT.0) Dg = 2 * Dg / Dqsc
       Dqsc =  2 * Dqsc / Dx**2
       Dqxt =  2 * Dqxt / Dx**2

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/pf_to_mr.F

@@ -107 +107 @@
                 mx_rain_cv(j,ibox,ilev)=0.
                 mx_snow_cv(j,ibox,ilev)=0.
-                if ((prec_frac(j,ibox,ilev) == 1.) .or.
-     &              (prec_frac(j,ibox,ilev) == 3.)) then
+                if ((prec_frac(j,ibox,ilev) .eq. 1.) .or.
+     &              (prec_frac(j,ibox,ilev) .eq. 3.)) then 
                     mx_rain_ls(j,ibox,ilev)=
      &                     (term4r_ls**(1./(1.+br/4.)))/rho
@@ -116 +116 @@
      &                     (term4g_ls**(1./(1.+bg/4.)))/rho
                 endif
-                if ((prec_frac(j,ibox,ilev) == 2.) .or.
-     &              (prec_frac(j,ibox,ilev) == 3.)) then
+                if ((prec_frac(j,ibox,ilev) .eq. 2.) .or.
+     &              (prec_frac(j,ibox,ilev) .eq. 3.)) then 
                     mx_rain_cv(j,ibox,ilev)=
      &                     (term4r_cv**(1./(1.+br/4.)))/rho

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/phys_cosp.F90

@@ -193 +193 @@
           cfg%Lrttov_sim,cfg%Lstats
 
-    if (overlaplmdz/=overlap) then
+    if (overlaplmdz.ne.overlap) then
        print*,'Attention overlaplmdz different de overlap lu dans namelist '
     endif
@@ -201 +201 @@
 
 !!! Ici on modifie les cles logiques pour les outputs selon les champs actives dans les .xml
-  if ((itap>1).and.(first_write))then
+  if ((itap.gt.1).and.(first_write))then
    
     IF (using_xios) call read_xiosfieldactive(cfg)
@@ -268 +268 @@
 
         do ip = 1, Npoints
-          if (fracTerLic(ip)>=0.5) then
+          if (fracTerLic(ip).ge.0.5) then
              gbx%land(ip) = 1.
           else

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/prec_scops.F

@@ -55 +55 @@
 
       cv_col = 0.05*ncol
-      if (cv_col == 0) cv_col=1
+      if (cv_col .eq. 0) cv_col=1
  
       do ilev=1,nlev
@@ -72 +72 @@
         flag_cv=0
         do ilev=1,nlev
-          if (frac_out(j,ibox,ilev) == 1) then
+          if (frac_out(j,ibox,ilev) .eq. 1) then 
             flag_ls=1
           endif
-          if (frac_out(j,ibox,ilev) == 2) then
+          if (frac_out(j,ibox,ilev) .eq. 2) then 
             flag_cv=1
           endif
         enddo !loop over nlev
-        if (flag_ls == 1) then
+        if (flag_ls .eq. 1) then
            frac_out_ls(j,ibox)=1
         endif
-        if (flag_cv == 1) then
+        if (flag_cv .eq. 1) then
            frac_out_cv(j,ibox)=1
         endif
@@ -93 +93 @@
         flag_cv=0
     
-        if (ls_p_rate(j,1) > 0.) then
+        if (ls_p_rate(j,1) .gt. 0.) then 
             do ibox=1,ncol ! possibility ONE
-                if (frac_out(j,ibox,1) == 1) then
+                if (frac_out(j,ibox,1) .eq. 1) then 
                     prec_frac(j,ibox,1) = 1
                     flag_ls=1
                 endif
             enddo ! loop over ncol
-            if (flag_ls == 0) then ! possibility THREE
+            if (flag_ls .eq. 0) then ! possibility THREE
                 do ibox=1,ncol
-                    if (frac_out(j,ibox,2) == 1) then
+                    if (frac_out(j,ibox,2) .eq. 1) then 
                         prec_frac(j,ibox,1) = 1
                         flag_ls=1
@@ -108 +108 @@
                 enddo ! loop over ncol
             endif
-        if (flag_ls == 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_ls(j,ibox) == 1) then
+        if (flag_ls .eq. 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_ls(j,ibox) .eq. 1) then 
             prec_frac(j,ibox,1) = 1
             flag_ls=1
@@ -116 +116 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls == 0) then ! possibility Five
+        if (flag_ls .eq. 0) then ! possibility Five
         do ibox=1,ncol
     !     prec_frac(j,1:ncol,1) = 1
@@ -125 +125 @@
        ! There is large scale precipitation
      
-        if (cv_p_rate(j,1) > 0.) then
+        if (cv_p_rate(j,1) .gt. 0.) then 
          do ibox=1,ncol ! possibility ONE
-          if (frac_out(j,ibox,1) == 2) then
-           if (prec_frac(j,ibox,1) == 0) then
+          if (frac_out(j,ibox,1) .eq. 2) then 
+           if (prec_frac(j,ibox,1) .eq. 0) then
         prec_frac(j,ibox,1) = 2
        else
@@ -136 +136 @@
       endif
         enddo ! loop over ncol
-        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
+        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
             prec_frac(j,ibox,1) = 2
             else
@@ -148 +148 @@
         enddo ! loop over ncol
         endif
-        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
+        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
             prec_frac(j,ibox,1) = 2
             else
@@ -160 +160 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv == 0) then  ! possibility Five
+        if (flag_cv .eq. 0) then  ! possibility Five
         do ibox=1,cv_col
-                if (prec_frac(j,ibox,1) == 0) then
+                if (prec_frac(j,ibox,1) .eq. 0) then
             prec_frac(j,ibox,1) = 2
             else
@@ -183 +183 @@
         flag_cv=0
     
-        if (ls_p_rate(j,ilev) > 0.) then
+        if (ls_p_rate(j,ilev) .gt. 0.) then 
          do ibox=1,ncol ! possibility ONE&TWO
-          if ((frac_out(j,ibox,ilev) == 1) .or.
-     &       ((prec_frac(j,ibox,ilev-1) == 1)
-     &       .or. (prec_frac(j,ibox,ilev-1) == 3))) then
+          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 
            prec_frac(j,ibox,ilev) = 1
            flag_ls=1
           endif
         enddo ! loop over ncol
-        if ((flag_ls == 0) .and. (ilev < nlev)) then ! possibility THREE
-        do ibox=1,ncol
-        if (frac_out(j,ibox,ilev+1) == 1) then
+        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 
             prec_frac(j,ibox,ilev) = 1
             flag_ls=1
@@ -200 +200 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls == 0) then ! possibility Four
-        do ibox=1,ncol
-        if (frac_out_ls(j,ibox) == 1) then
+        if (flag_ls .eq. 0) then ! possibility Four
+        do ibox=1,ncol
+        if (frac_out_ls(j,ibox) .eq. 1) then 
             prec_frac(j,ibox,ilev) = 1
             flag_ls=1
@@ -208 +208 @@
         enddo ! loop over ncol
         endif
-        if (flag_ls == 0) then ! possibility Five
+        if (flag_ls .eq. 0) then ! possibility Five
         do ibox=1,ncol
 !     prec_frac(j,1:ncol,ilev) = 1
@@ -216 +216 @@
       endif ! There is large scale precipitation
     
-        if (cv_p_rate(j,ilev) > 0.) then
+        if (cv_p_rate(j,ilev) .gt. 0.) then 
          do ibox=1,ncol ! possibility ONE&TWO
-          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
+          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
          prec_frac(j,ibox,ilev) = 2
         else
@@ -229 +229 @@
         endif
        enddo ! loop over ncol
-        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
+        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
             prec_frac(j,ibox,ilev) = 2
             else
@@ -241 +241 @@
         enddo ! loop over ncol
         endif
-        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
+        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
             prec_frac(j,ibox,ilev) = 2
             else
@@ -253 +253 @@
         enddo ! loop over ncol
         endif
-        if (flag_cv == 0) then  ! possibility Five
+        if (flag_cv .eq. 0) then  ! possibility Five 
         do ibox=1,cv_col
-                if (prec_frac(j,ibox,ilev) == 0) then
+                if (prec_frac(j,ibox,ilev) .eq. 0) then
             prec_frac(j,ibox,ilev) = 2
             else

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/predict_mom07.F90

@@ -7 +7 @@
             implicit none
         
-            real(kind=8) :: a1,a2,a3,b1,b2,b3,c1,c2,c3
-            real(kind=8) :: m2,tc,n,m,a_,b_,c_,A,B,C,n2
+            real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3
+            real*8 :: m2,tc,n,m,a_,b_,c_,A,B,C,n2
         
             a1=      13.6078
@@ -30 +30 @@
             
         ! predict m from m2 and tc
-                if(m2/=-9999) then
+                if(m2.ne.-9999) then 
                 m=A*exp(B*tc)*m2**C
                 endif
         ! get m2 if mass-dimension relationship not proportional to D**2
-                if(m2==-9999) then
+                if(m2.eq.-9999) then 
                 m2=(m/(A*exp(B*tc)))**(1.0/C)    
                 endif

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator.F90

@@ -91 +91 @@
 
   real undef
-  real(kind=8), dimension(nprof,ngate), intent(in) :: &
+  real*8, dimension(nprof,ngate), intent(in) :: &
     hgt_matrix, p_matrix,t_matrix,rh_matrix
 
-  real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
-  real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
-  real(kind=8), dimension(hp%nhclass,nprof,ngate), intent(inout)    :: Np_matrix
+  real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
+  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
+  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout)    :: Np_matrix
 
 ! ----- OUTPUTS -----
-  real(kind=8), dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
+  real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
        g_to_vol,dBZe,a_to_vol
 
 ! ----- OPTIONAL -----
-  real(kind=8), optional, dimension(nprof,ngate) :: &
+  real*8, optional, dimension(nprof,ngate) :: &
   g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
                            ! the same gaseous absorption in different calls. Optional to allow compatibility
@@ -112 +112 @@
 
   real, parameter :: one_third = 1.0/3.0
-  real(kind=8) :: t_kelvin
+  real*8 :: t_kelvin
   integer :: &
   phase, & ! 0=liquid, 1=ice
@@ -118 +118 @@
 
   logical :: hydro      ! true=hydrometeor in vol, false=none
-  real(kind=8) :: &
+  real*8 :: &
   rho_a, &   ! air density (kg m^-3)
   gases      ! function: 2-way gas atten (dB/km)
 
-  real(kind=8), dimension(:), allocatable :: &
+  real*8, dimension(:), allocatable :: &
   Di, Deq, &   ! discrete drop sizes (um)
   Ni, &        ! discrete concentrations (cm^-3 um^-1)
   rhoi         ! discrete densities (kg m^-3)
 
-  real(kind=8), dimension(nprof, ngate) :: &
+  real*8, dimension(nprof, ngate) :: &
   z_vol, &      ! effective reflectivity factor (mm^6/m^3)
   z_ray, &                      ! reflectivity factor, Rayleigh only (mm^6/m^3)
@@ -135 +135 @@
 
   integer,parameter :: KR8 = selected_real_kind(15,300)
-  real(kind=8), parameter :: xx = -1.0_KR8
-  real(kind=8),  dimension(:), allocatable :: xxa
-  real(kind=8) :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
-  real(kind=8) :: half_a_atten_current,half_a_atten_above
-  real(kind=8) :: half_g_atten_current,half_g_atten_above
-  integer(kind=4) :: tp, i, j, k, pr, itt, iff
-
-  real(kind=8)    step,base, Np
-  integer(kind=4) iRe_type,n,max_bin
+  real*8, parameter :: xx = -1.0_KR8
+  real*8,  dimension(:), allocatable :: xxa
+  real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
+  real*8 :: half_a_atten_current,half_a_atten_above
+  real*8 :: half_g_atten_current,half_g_atten_above
+  integer*4 :: tp, i, j, k, pr, itt, iff
+
+  real*8    step,base, Np
+  integer*4 iRe_type,n,max_bin
 
   integer   start_gate,end_gate,d_gate
@@ -207 +207 @@
             itt = infind(hp%mt_tti,t_kelvin)
           endif
-          if (re_matrix(tp,pr,k)==0) then
+          if (re_matrix(tp,pr,k).eq.0) then
             call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, &
               hp%dtype(tp),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
@@ -221 +221 @@
 
           iRe_type=1
-          if(Re>0) then
+          if(Re.gt.0) then
             ! determine index in to scale LUT
             !
@@ -232 +232 @@
             base=hp%base_list(n+1)
             iRe_type=Re/step
-            if (iRe_type<1) iRe_type=1
+            if (iRe_type.lt.1) iRe_type=1
 
             Re=step*(iRe_type+0.5)      ! set value of Re to closest value allowed in LUT.
@@ -238 +238 @@
 
             ! make sure iRe_type is within bounds
-            if (iRe_type>=nRe_types) then
+            if (iRe_type.ge.nRe_types) then
 !               write(*,*) 'Warning: size of Re exceed value permitted ', &
 !                    'in Look-Up Table (LUT).  Will calculate. '

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_init.F90

@@ -74 +74 @@
 ! ----- INTERNAL -----  
   integer :: i,j
-  real(kind=8)  :: delt, deltp
+  real*8  :: delt, deltp
         
     !

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_types.F90

@@ -12 +12 @@
   integer, parameter ::       &
   nd = 85               ! number of discrete particles used in construction DSDs
-  real(kind=8), parameter ::        &
+  real*8, parameter ::        &
   dmin = 0.1                 ,& ! min size of discrete particle
   dmax = 10000.                 ! max size of discrete particle
@@ -36 +36 @@
   
     ! variables used to store hydrometeor "default" properties
-    real(kind=8),  dimension(maxhclass) :: p1,p2,p3,dmin,dmax,apm,bpm,rho
+    real*8,  dimension(maxhclass) :: p1,p2,p3,dmin,dmax,apm,bpm,rho
     integer, dimension(maxhclass) :: dtype,col,cp,phase
   
     ! Radar properties
-    real(kind=8)  :: freq,k2
+    real*8  :: freq,k2
     integer :: nhclass      ! number of hydrometeor classes in use
     integer :: use_gas_abs, do_ray
@@ -56 +56 @@
     logical, dimension(maxhclass,nRe_types) :: N_scale_flag
     logical, dimension(maxhclass,mt_ntt,nRe_types) :: Z_scale_flag,Z_scale_added_flag
-    real(kind=8),  dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled
-    real(kind=8),  dimension(maxhclass,nd,nRe_types) :: fc, rho_eff
+    real*8,  dimension(maxhclass,mt_ntt,nRe_types) :: Ze_scaled,Zr_scaled,kr_scaled
+    real*8,  dimension(maxhclass,nd,nRe_types) :: fc, rho_eff
 
     ! used to determine Re index
-    real(kind=8)  :: step_list(Re_MAX_BIN),base_list(Re_MAX_BIN)
+    real*8  :: step_list(Re_MAX_BIN),base_list(Re_MAX_BIN)
   
     ! used to determine temperature index
-    real(kind=8) :: &
+    real*8 :: &
         mt_ttl(cnt_liq), &  ! liquid temperatures (K)
         mt_tti(cnt_ice)     ! ice temperatures (K)
 
-    real(kind=8) :: D(nd) ! set of discrete diameters used to represent DSDs
+    real*8 :: D(nd) ! set of discrete diameters used to represent DSDs
 
   end type class_param

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/scops.F

@@ -133 +133 @@
       enddo
 
-      if (ncolprint/=0) then
+      if (ncolprint.ne.0) then
         write (6,'(a)') 'frac_out_pp_rev:'
           do j=1,npoints,1000
@@ -145 +145 @@
         write (6,'(I3)') ncol
       endif
-      if (ncolprint/=0) then
+      if (ncolprint.ne.0) then
         write (6,'(a)') 'last_frac_pp:'
           do j=1,npoints,1000
@@ -161 +161 @@
       
       !loop over vertical levels
-      DO ilev = 1,nlev
+      DO 200 ilev = 1,nlev
                                   
 !     Initialise threshold
 
-        IF (ilev==1) then
+        IF (ilev.eq.1) then
           ! If max overlap 
-          IF (overlap==1) then
+          IF (overlap.eq.1) then
             ! select pixels spread evenly
             ! across the gridbox
@@ -187 +187 @@
               enddo
             ENDIF
-            IF (ncolprint/=0) then
+            IF (ncolprint.ne.0) then
               write (6,'(a)') 'threshold_nsf2:'
                 do j=1,npoints,1000
@@ -197 +197 @@
         ENDIF
 
-        IF (ncolprint/=0) then
+        IF (ncolprint.ne.0) then
             write (6,'(a)') 'ilev:'
             write (6,'(I2)') ilev
@@ -206 +206 @@
           ! All versions
           do j=1,npoints
-            if (boxpos(j,ibox)<=conv(j,ilev)) then
+            if (boxpos(j,ibox).le.conv(j,ilev)) then
               maxocc(j,ibox) = 1
             else
@@ -214 +214 @@
 
           ! Max overlap
-          if (overlap==1) then
+          if (overlap.eq.1) then 
             do j=1,npoints
               threshold_min(j,ibox)=conv(j,ilev)
@@ -222 +222 @@
 
           ! Random overlap
-          if (overlap==2) then
+          if (overlap.eq.2) then 
             do j=1,npoints
               threshold_min(j,ibox)=conv(j,ilev)
@@ -230 +230 @@
 
           ! Max/Random overlap
-          if (overlap==3) then
+          if (overlap.eq.3) then 
             do j=1,npoints
               threshold_min(j,ibox)=max(conv(j,ilev),
      &          min(tca(j,ilev-1),tca(j,ilev)))
               if (threshold(j,ibox)
-     &          <min(tca(j,ilev-1),tca(j,ilev))
-     &          .and.(threshold(j,ibox)>conv(j,ilev))) then
+     &          .lt.min(tca(j,ilev-1),tca(j,ilev))
+     &          .and.(threshold(j,ibox).gt.conv(j,ilev))) then
                    maxosc(j,ibox)= 1
               else
@@ -276 +276 @@
            DO ibox=1,ncol
              do j=1,npoints 
-               if (tca(j,ilev)>threshold(j,ibox)) then
+               if (tca(j,ilev).gt.threshold(j,ibox)) then
                frac_out(j,ibox,ilev)=1
                else
@@ -289 +289 @@
            DO ibox=1,ncol
              do j=1,npoints 
-                if (threshold(j,ibox)<=conv(j,ilev)) then
+                if (threshold(j,ibox).le.conv(j,ilev)) then
                     ! = 2 IF threshold le conv(j)
                     frac_out(j,ibox,ilev) = 2 
@@ -302 +302 @@
 !         from last level next time round
 
-          if (ncolprint/=0) then
+          if (ncolprint.ne.0) then
 
             do j=1,npoints ,1000
@@ -331 +331 @@
           endif
 
-      END DO    !loop over nlev
+200   CONTINUE    !loop over nlev
 
 

LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/zeff.F90

@@ -35 +35 @@
   integer, intent(in) :: ice, xr
   integer, intent(in) :: nsizes
-  real(kind=8), intent(in) :: freq,D(nsizes),N(nsizes),tt,qe(nsizes), &
+  real*8, intent(in) :: freq,D(nsizes),N(nsizes),tt,qe(nsizes), &
     qs(nsizes), rho_e(nsizes)
-  real(kind=8), intent(inout) :: k2
+  real*8, intent(inout) :: k2
   
 ! ----- OUTPUTS -----
-  real(kind=8), intent(out) :: z_eff,z_ray,kr
+  real*8, intent(out) :: z_eff,z_ray,kr
     
 ! ----- INTERNAL -----
   integer :: &
   correct_for_rho        ! correct for density flag
-  real(kind=8), dimension(nsizes) :: &
+  real*8, dimension(nsizes) :: &
   D0, &                  ! D in (m)
   N0, &                  ! N in m^-3 m^-1
@@ -53 +53 @@
   rho_ice, &             ! bulk density ice (kg m^-3)
   f                 ! ice fraction
-  real(kind=8), dimension(nsizes) :: xtemp
-  real(kind=8) :: &
+  real*8, dimension(nsizes) :: xtemp
+  real*8 :: &
   wl, &                  ! wavelength (m)
   cr                            ! kr(dB/km) = cr * kr(1/km)
-  complex(kind=8) :: &
+  complex*16 :: &
   m                 ! complex index of refraction of bulk form
-  complex(kind=8), dimension(nsizes) :: &
+  complex*16, dimension(nsizes) :: &
   m0                ! complex index of refraction
   
-  integer(kind=4) :: i,one
-  real(kind=8) :: pi
-  real(kind=8) :: eta_sum, eta_mie, const, z0_eff, z0_ray, k_sum, &
+  integer*4 :: i,one
+  real*8 :: pi
+  real*8 :: eta_sum, eta_mie, const, z0_eff, z0_ray, k_sum, &
             n_r, n_i, dqv(1), dqsc, dg, dph(1)
-  integer(kind=4) :: err
-  complex(kind=8) :: Xs1(1), Xs2(1)
+  integer*4 :: err
+  complex*16 :: Xs1(1), Xs2(1)
 
   one=1
@@ -113 +113 @@
       call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), &
         dg, xs1, xs2, dph, err)
-    END DO
+    end do
     
   else