Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver.F90	(revision 5158)
@@ -311,5 +311,5 @@
   tstart = omp_get_wtime() 
 #endif
-  do jrepeat = 1,driver_config%nrepeat
+  DO jrepeat = 1,driver_config%nrepeat
     
     if (driver_config%do_parallel) then
@@ -321,5 +321,5 @@
      
       !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME)
-      do jblock = 1, nblock
+      DO jblock = 1, nblock
         ! Specify the range of columns to process.
         istartcol = (jblock-1) * driver_config%nblocksize &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_config.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_config.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_config.F90	(revision 5158)
@@ -400,5 +400,5 @@
     ! wavelength bound, noting that the number of diagnostics is one
     ! fewer than the number of valid bounds
-    do jdiag = 0,NMaxSpectralDiag
+    DO jdiag = 0,NMaxSpectralDiag
       if (this%sw_diag_wavelength_bound(jdiag+1) < 0.0_jprb) then
         this%n_sw_diag = max(0,jdiag-1)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_read_input.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_read_input.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_read_input.F90	(revision 5158)
@@ -563,5 +563,5 @@
 
     ! Loop through all radiatively important gases
-    do jgas = 1,NMaxGases
+    DO jgas = 1,NMaxGases
       if (jgas == IH2O) then
         if (file%exists('q')) then
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver.F90	(revision 5158)
@@ -369,5 +369,5 @@
   tstart = omp_get_wtime()
 #endif
-  do jrepeat = 1,driver_config%nrepeat
+  DO jrepeat = 1,driver_config%nrepeat
 
 !    if (driver_config%do_parallel) then
@@ -379,5 +379,5 @@
 
       !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME)
-      do jblock = 1, nblock
+      DO jblock = 1, nblock
         ! Specify the range of columns to process.
         istartcol = (jblock-1) * driver_config%nblocksize &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver_blocked.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver_blocked.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver_blocked.F90	(revision 5158)
@@ -371,5 +371,5 @@
   tstart = omp_get_wtime()
 #endif
-  do jrepeat = 1,driver_config%nrepeat
+  DO jrepeat = 1,driver_config%nrepeat
 
 !    if (driver_config%do_parallel) then
@@ -378,5 +378,5 @@
       !$OMP PARALLEL DO SCHEDULE(DYNAMIC,1)&
       !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB)
-      do jrl=1,ncol,nproma
+      DO jrl=1,ncol,nproma
         ibeg=jrl
         iend=min(ibeg+nproma-1,ncol)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ifs_blocking.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ifs_blocking.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ifs_blocking.F90	(revision 5158)
@@ -329,6 +329,6 @@
   !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
   !$OMP&PRIVATE(IB,IFLD)
-  do ib=1,ngpblks
-    do ifld=1,ifs_config%ifldstot
+  DO ib=1,ngpblks
+    DO ifld=1,ifs_config%ifldstot
       zrgp(:,ifld,ib) = 0._jprb
     enddo
@@ -351,5 +351,5 @@
     !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
     !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JAER,JOFF,JLEV,JALB)
-    do jrl=1,ncol,nproma
+    DO jrl=1,ncol,nproma
 
       ibeg=jrl
@@ -363,5 +363,5 @@
       zrgp(1:il,ifs_config%iamu0,ib)  =  single_level%cos_sza(ibeg:iend)   ! cosine of solar zenith ang (mu0)
 
-      do jemiss=1,yderad%nlwemiss
+      DO jemiss=1,yderad%nlwemiss
         zrgp(1:il,ifs_config%iemiss+jemiss-1,ib)  =  single_level%lw_emissivity(ibeg:iend,jemiss)
       enddo
@@ -378,19 +378,19 @@
       ! zrgp(1:il,islon,ib)    = ???
 
-      do jalb=1,yderad%nsw
+      DO jalb=1,yderad%nsw
         zrgp(1:il,ifs_config%iald+jalb-1,ib)  =  single_level%sw_albedo(ibeg:iend,jalb)
       enddo
 
       if (allocated(single_level%sw_albedo_direct)) then
-        do jalb=1,yderad%nsw
+        DO jalb=1,yderad%nsw
           zrgp(1:il,ifs_config%ialp+jalb-1,ib)  =  single_level%sw_albedo_direct(ibeg:iend,jalb)
         end do
       else
-        do jalb=1,yderad%nsw
+        DO jalb=1,yderad%nsw
           zrgp(1:il,ifs_config%ialp+jalb-1,ib)  =  single_level%sw_albedo(ibeg:iend,jalb)
         end do
       end if
       
-      do jlev=1,nlev
+      DO jlev=1,nlev
         zrgp(1:il,ifs_config%iti+jlev-1,ib)   = temperature_fl(ibeg:iend,jlev) ! full level temperature
         zrgp(1:il,ifs_config%ipr+jlev-1,ib)   = pressure_fl(ibeg:iend,jlev) ! full level pressure
@@ -398,5 +398,5 @@
       enddo
 
-      do jlev=1,nlev
+      DO jlev=1,nlev
         zrgp(1:il,ifs_config%iwv+jlev-1,ib)   = gas%mixing_ratio(ibeg:iend,jlev,IH2O) ! this is already in MassMixingRatio units
         if (rad_config%do_clouds) then
@@ -421,6 +421,6 @@
       if (yderad%naermacc == 1) then
         joff=ifs_config%iaero
-        do jaer=1,rad_config%n_aerosol_types
-          do jlev=1,nlev
+        DO jaer=1,rad_config%n_aerosol_types
+          DO jlev=1,nlev
             zrgp(1:il,joff,ib) = aerosol%mixing_ratio(ibeg:iend,jlev,jaer)
             joff=joff+1
@@ -429,5 +429,5 @@
       endif
 
-      do jlev=1,nlev+1
+      DO jlev=1,nlev+1
         ! zrgp(1:il,ihpr+jlev-1,ib)  = ???
         zrgp(1:il,ifs_config%iaprs+jlev-1,ib) = thermodynamics%pressure_hl(ibeg:iend,jlev)
@@ -453,10 +453,10 @@
       ! local workaround variables for standalone input files
       if (rad_config%do_clouds) then
-        do jlev=1,nlev
+        DO jlev=1,nlev
           ! missing full-level temperature and pressure as well as land-sea-mask
           zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = cloud%re_liq(ibeg:iend,jlev)
           zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = cloud%re_ice(ibeg:iend,jlev)
         enddo
-        do jlev=1,nlev-1
+        DO jlev=1,nlev-1
           ! for the love of it, I can't figure this one out. Probably to do with
           ! my crude approach of setting PGEMU?
@@ -465,10 +465,10 @@
         if(present(iseed)) iseed(1:il,ib) = single_level%iseed(ibeg:iend)
       else
-        do jlev=1,nlev
+        DO jlev=1,nlev
           ! missing full-level temperature and pressure as well as land-sea-mask
           zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = 0._jprb
           zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = 0._jprb
         enddo
-        do jlev=1,nlev-1
+        DO jlev=1,nlev-1
           zrgp(1:il,ifs_config%ioverlap+jlev-1,ib) = 0._jprb
         enddo
@@ -531,5 +531,5 @@
     !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
     !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JLEV,JG)
-    do jrl=1,ncol,nproma
+    DO jrl=1,ncol,nproma
       ibeg=jrl
       iend=min(ibeg+nproma-1,ncol)
@@ -537,5 +537,5 @@
       ib=(jrl-1)/nproma+1
 
-      do jlev=1,nlev+1
+      DO jlev=1,nlev+1
         flux%sw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrso+jlev-1,ib)
         flux%lw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrth+jlev-1,ib)
@@ -561,5 +561,5 @@
       emissivity_out(ibeg:iend) = zrgp(1:il,ifs_config%iemit,ib)
       if (yradiation%yrerad%lapproxswupdate) then
-        do jg=1,yradiation%yrerad%nsw
+        DO jg=1,yradiation%yrerad%nsw
           flux_diffuse_band(ibeg:iend,jg) = zrgp(1:il,ifs_config%iswdiffuseband+jg-1,ib)
           flux_direct_band(ibeg:iend,jg) = zrgp(1:il,ifs_config%iswdirectband+jg-1,ib)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_cloud_generator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_cloud_generator.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_cloud_generator.F90	(revision 5158)
@@ -57,6 +57,6 @@
        &  od_scaling, total_cloud_cover)
 
-  do jlev = 1,nlev
-    do jcol = 1,ncol
+  DO jlev = 1,nlev
+    DO jcol = 1,ncol
       write(6,'(f5.2,a)',advance='no') od_scaling(jcol,jlev), ' '
     end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_random_number_generator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_random_number_generator.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_random_number_generator.F90	(revision 5158)
@@ -29,5 +29,5 @@
   call random_number_generator%initialize(itype=IRngType, iseed=212075152, &
        &                                  nmaxstreams=streammax)
-  do jl = 1,1
+  DO jl = 1,1
     print *, 'initial_state = [ ', int(random_number_generator%istate(1:streammax),jpib), ' ]'
     call random_number_generator%uniform_distribution(vec)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_spartacus_math.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_spartacus_math.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_spartacus_math.F90	(revision 5158)
@@ -45,10 +45,10 @@
 
   write(*,*) 'A ='
-  do j = 1,m
+  DO j = 1,m
      write(*,*) A(1,j,:)
   end do
 
   write(*,*) 'B ='
-  do j = 1,m
+  DO j = 1,m
      write(*,*) B(1,j,:)
   end do
@@ -60,5 +60,5 @@
 
   write(*,*) 'C = A*B ='
-  do j = 1,m
+  DO j = 1,m
      write(*,*) C(1,j,:)
   end do
@@ -70,5 +70,5 @@
 
   write(*,*) 'C = A\B ='
-  do j = 1,m
+  DO j = 1,m
      write(*,*) C(1,j,:)
   end do
@@ -78,5 +78,5 @@
 
   write(*,*) 'expm(A) ='
-  do j = 1,m
+  DO j = 1,m
      write(*,*) A(1,j,:)
   end do
@@ -93,5 +93,5 @@
     
     write(*,*) 'fast_expm(A) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
@@ -106,5 +106,5 @@
     
     write(*,*) 'expm(A) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
@@ -115,5 +115,5 @@
     
     write(*,*) 'expm(zeros) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
@@ -124,5 +124,5 @@
     
     write(*,*) 'fast_expm(A) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
@@ -140,5 +140,5 @@
     
     write(*,*) 'expm(A) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
@@ -149,5 +149,5 @@
     
     write(*,*) 'expm(zeros) = '
-    do j = 1,m
+    DO j = 1,m
       write(*,*) A(1,j,:)
     end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/fcttre.func.h
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/fcttre.func.h	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/fcttre.func.h	(revision 5158)
@@ -170,7 +170,7 @@
 FOEELIQ2ICE( PTARE )  = EXP(      (PTARE-RTT)*(R3LES/(PTARE-R4LES) - R3IES/(PTARE-R4IES)))
 FOELSON( PTARE ) = EXP( -6096.9385_JPRB/PTARE + 21.2409642_JPRB &
-	             - 2.711193E-2_JPRB * PTARE    &
+                 - 2.711193E-2_JPRB * PTARE    &
                      + 1.673952E-5_JPRB * PTARE**2 &
-		     + 2.433502_JPRB * LOG(PTARE))
+             + 2.433502_JPRB * LOG(PTARE))
 
 REAL(KIND=JPRB) :: FOEEWM_V,FOEEWMCU_V,FOELES_V,FOEIES_V
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90	(revision 5158)
@@ -54,8 +54,8 @@
 
       ! Loop over level
-     do jlev = istartlev,iendlev
+     DO jlev = istartlev,iendlev
 
         ! Loop over column
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
 
           ! Compute relative humidity with respect to liquid
@@ -71,5 +71,5 @@
           od_aerosol_mono = 0.0_jprb
 
-          do jtype = 1,config%n_aerosol_types
+          DO jtype = 1,config%n_aerosol_types
             ! Add the optical depth for this aerosol type to the total for all aerosols.  
             ! Note that the following expressions are array-wise
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_lw.F90	(revision 5158)
@@ -83,5 +83,5 @@
     ! also the "source", which is the upwelling flux due to emission
     ! below that level
-    do jlev = nlev,1,-1
+    DO jlev = nlev,1,-1
       ! Next loop over columns. We could do this by indexing the
       ! entire inner dimension as follows, e.g. for the first line:
@@ -90,5 +90,5 @@
       ! routine by a factor of 2!  Rather, we do it with an explicit
       ! loop.
-      do jcol = 1,ncol
+      DO jcol = 1,ncol
         ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
         inv_denominator(jcol,jlev) = 1.0_jprb &
@@ -114,6 +114,6 @@
     ! Work back down through the atmosphere computing the fluxes at
     ! each half-level
-    do jlev = 1,nlev
-      do jcol = 1,ncol
+    DO jlev = 1,nlev
+      DO jcol = 1,ncol
         ! Shonk & Hogan (2008) Eq 14 (after simplification):
         flux_dn(jcol,jlev+1) &
@@ -201,8 +201,8 @@
     ! also the "source", which is the upwelling flux due to emission
     ! below that level
-    do jlev = nlev,i_cloud_top,-1
+    DO jlev = nlev,i_cloud_top,-1
       if (is_clear_sky_layer(jlev)) then
         ! Reflectance of this layer is zero, simplifying the expression
-        do jcol = 1,ncol
+        DO jcol = 1,ncol
           albedo(jcol,jlev) = transmittance(jcol,jlev)*transmittance(jcol,jlev)*albedo(jcol,jlev+1)
           source(jcol,jlev) = source_up(jcol,jlev) &
@@ -212,5 +212,5 @@
       else
         ! Loop over columns; explicit loop seems to be faster
-        do jcol = 1,ncol
+        DO jcol = 1,ncol
           ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
           inv_denominator(jcol,jlev) = 1.0_jprb &
@@ -231,5 +231,5 @@
     flux_up(:,i_cloud_top) = source(:,i_cloud_top) &
          &                 + albedo(:,i_cloud_top)*flux_dn(:,i_cloud_top)
-    do jlev = i_cloud_top-1,1,-1
+    DO jlev = i_cloud_top-1,1,-1
       flux_up(:,jlev) = transmittance(:,jlev)*flux_up(:,jlev+1) + source_up(:,jlev)
     end do
@@ -237,7 +237,7 @@
     ! Work back down through the atmosphere from cloud top computing
     ! the fluxes at each half-level
-    do jlev = i_cloud_top,nlev
+    DO jlev = i_cloud_top,nlev
       if (is_clear_sky_layer(jlev)) then
-        do jcol = 1,ncol
+        DO jcol = 1,ncol
           flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) &
                &               + source_dn(jcol,jlev)
@@ -246,5 +246,5 @@
         end do
       else
-        do jcol = 1,ncol
+        DO jcol = 1,ncol
           ! Shonk & Hogan (2008) Eq 14 (after simplification):
           flux_dn(jcol,jlev+1) &
@@ -309,6 +309,6 @@
 ! Added for DWD (2020)
 !NEC$ outerloop_unroll(8)
-    do jlev = 1,nlev
-      do jcol = 1,ncol
+    DO jlev = 1,nlev
+      DO jcol = 1,ncol
         flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) + source_dn(jcol,jlev)
       end do
@@ -322,6 +322,6 @@
 ! Added for DWD (2020)
 !NEC$ outerloop_unroll(8)
-    do jlev = nlev,1,-1
-      do jcol = 1,ncol
+    DO jlev = nlev,1,-1
+      DO jcol = 1,ncol
         flux_up(jcol,jlev) = transmittance(jcol,jlev)*flux_up(jcol,jlev+1) + source_up(jcol,jlev)
       end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_sw.F90	(revision 5158)
@@ -84,5 +84,5 @@
     ! half-level by working down through the atmosphere
     flux_dn_direct(:,1) = incoming_toa
-    do jlev = 1,nlev
+    DO jlev = 1,nlev
       flux_dn_direct(:,jlev+1) = flux_dn_direct(:,jlev)*trans_dir_dir(:,jlev)
     end do
@@ -100,5 +100,5 @@
 ! Added for DWD (2020)
 !NEC$ outerloop_unroll(8)
-    do jlev = nlev,1,-1
+    DO jlev = nlev,1,-1
       ! Next loop over columns. We could do this by indexing the
       ! entire inner dimension as follows, e.g. for the first line:
@@ -107,5 +107,5 @@
       ! routine by a factor of 2!  Rather, we do it with an explicit
       ! loop.
-      do jcol = 1,ncol
+      DO jcol = 1,ncol
         ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
         inv_denominator(jcol,jlev) = 1.0_jprb / (1.0_jprb-albedo(jcol,jlev+1)*reflectance(jcol,jlev))
@@ -132,6 +132,6 @@
 ! Added for DWD (2020)
 !NEC$ outerloop_unroll(8)
-    do jlev = 1,nlev
-      do jcol = 1,ncol
+    DO jlev = 1,nlev
+      DO jcol = 1,ncol
         ! Shonk & Hogan (2008) Eq 14 (after simplification):
         flux_dn_diffuse(jcol,jlev+1) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90	(revision 5158)
@@ -243,5 +243,5 @@
 
         if (ao%use_hydrophilic) then
-          do jtype = 1,n_type_philic
+          DO jtype = 1,n_type_philic
             ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype))
             ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) &
@@ -266,5 +266,5 @@
 
         if (ao%use_hydrophilic) then
-          do jtype = 1,n_type_philic
+          DO jtype = 1,n_type_philic
             ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype))
             ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) &
@@ -280,5 +280,5 @@
       if (allocated(ao%wavelength_mono)) then
         ! Monochromatic wavelengths also required
-        do jwl = 1,nmono
+        DO jwl = 1,nmono
           ! Wavelength (m) to wavenumber (cm-1)
           wavenumber_target = 0.01_jprb / ao%wavelength_mono(jwl)
@@ -292,5 +292,5 @@
           else
             iwn = 1
-            do while (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1)
+            DO while (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1)
               iwn = iwn + 1
             end do
@@ -419,5 +419,5 @@
 
       if (ao%use_hydrophilic) then
-        do jtype = 1,ao%n_type_philic
+        DO jtype = 1,ao%n_type_philic
           ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype))
           ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype) &
@@ -451,5 +451,5 @@
 
       if (ao%use_hydrophilic) then
-        do jtype = 1,ao%n_type_philic
+        DO jtype = 1,ao%n_type_philic
           ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype))
           ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype) &
@@ -579,5 +579,5 @@
       ! Aerosol mixing ratios have been provided
 
-      do jtype = 1,config%n_aerosol_types
+      DO jtype = 1,config%n_aerosol_types
         if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
           write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
@@ -612,5 +612,5 @@
 
       ! Loop over column
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
 
         ! Reset temporary arrays
@@ -622,5 +622,5 @@
         scat_g_lw_aerosol = 0.0_jprb
 
-        do jlev = istartlev,iendlev
+        DO jlev = istartlev,iendlev
           ! Compute relative humidity with respect to liquid
           ! saturation and the index to the relative-humidity index of
@@ -634,5 +634,5 @@
         end do
 
-        do jtype = 1,config%n_aerosol_types
+        DO jtype = 1,config%n_aerosol_types
           itype = ao%itype(jtype)
 
@@ -643,7 +643,7 @@
           ! dimension being spectral band.
           if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
-            do jlev = istartlev,iendlev
+            DO jlev = istartlev,iendlev
               mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype)
-              do jband = 1,config%n_bands_sw
+              DO jband = 1,config%n_bands_sw
                 local_od_sw = factor(jlev) * mixing_ratio &
                      &  * ao%mass_ext_sw_phobic(jband,itype)
@@ -656,5 +656,5 @@
               end do
               if (config%do_lw_aerosol_scattering) then
-                do jband = 1,config%n_bands_lw
+                DO jband = 1,config%n_bands_lw
                   local_od_lw = factor(jlev) * mixing_ratio &
                        &  * ao%mass_ext_lw_phobic(jband,itype)
@@ -670,5 +670,5 @@
                 ! weight the optical depth by the single scattering
                 ! co-albedo
-                do jband = 1,config%n_bands_lw
+                DO jband = 1,config%n_bands_lw
                   od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) &
                        &  + factor(jlev) * mixing_ratio &
@@ -682,8 +682,8 @@
             ! Hydrophilic aerosols require the look-up tables to
             ! be indexed with irh
-            do jlev = istartlev,iendlev
+            DO jlev = istartlev,iendlev
               mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype)
               irh = irhs(jlev)
-              do jband = 1,config%n_bands_sw
+              DO jband = 1,config%n_bands_sw
                 local_od_sw = factor(jlev) * mixing_ratio &
                      &  * ao%mass_ext_sw_philic(jband,irh,itype)
@@ -696,5 +696,5 @@
               end do
               if (config%do_lw_aerosol_scattering) then
-                do jband = 1,config%n_bands_lw
+                DO jband = 1,config%n_bands_lw
                   local_od_lw = factor(jlev) * mixing_ratio &
                        &  * ao%mass_ext_lw_philic(jband,irh,itype)
@@ -710,5 +710,5 @@
                 ! weight the optical depth by the single scattering
                 ! co-albedo
-                do jband = 1,config%n_bands_lw
+                DO jband = 1,config%n_bands_lw
                   od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) &
                        &  + factor(jlev) * mixing_ratio &
@@ -740,6 +740,6 @@
 
           ! We can assume the band and g-point indices are the same
-          do jlev = 1,nlev
-            do jg = 1,config%n_g_sw
+          DO jlev = 1,nlev
+            DO jg = 1,config%n_g_sw
               local_scat = ssa_sw(jg,jlev,jcol)*od_sw(jg,jlev,jcol) + scat_sw_aerosol(jg,jlev)
               od_sw(jg,jlev,jcol) = od_sw(jg,jlev,jcol) + od_sw_aerosol(jg,jlev)
@@ -751,6 +751,6 @@
         else
 
-          do jlev = 1,nlev
-            do jg = 1,config%n_g_sw
+          DO jlev = 1,nlev
+            DO jg = 1,config%n_g_sw
               ! Need to map between bands and g-points
               iband = config%i_band_from_reordered_g_sw(jg)
@@ -781,6 +781,6 @@
                &                             scat_lw_aerosol, scat_g_lw_aerosol)
 
-          do jlev = istartlev,iendlev
-            do jg = 1,config%n_g_lw
+          DO jlev = istartlev,iendlev
+            DO jg = 1,config%n_g_lw
               iband = config%i_band_from_reordered_g_lw(jg)
               local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband,jlev)
@@ -802,12 +802,12 @@
           if (config%do_cloud_aerosol_per_lw_g_point) then
             ! We can assume band and g-point indices are the same
-            do jlev = istartlev,iendlev
-              do jg = 1,config%n_g_lw
+            DO jlev = istartlev,iendlev
+              DO jg = 1,config%n_g_lw
                 od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) + od_lw_aerosol(jg,jlev)
               end do
             end do
           else
-            do jlev = istartlev,iendlev
-              do jg = 1,config%n_g_lw
+            DO jlev = istartlev,iendlev
+              DO jg = 1,config%n_g_lw
                 od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
                      &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev)
@@ -901,9 +901,9 @@
 
       ! Loop over position
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
 ! Added for DWD (2020)
 !NEC$ forced_collapse
-        do jlev = istartlev,iendlev
-          do jb = 1,config%n_bands_sw
+        DO jlev = istartlev,iendlev
+          DO jb = 1,config%n_bands_sw
             od_sw_aerosol(jb,jlev) = aerosol%od_sw(jb,jlev,jcol)
             scat_sw_aerosol(jb,jlev) = aerosol%ssa_sw(jb,jlev,jcol) * od_sw_aerosol(jb,jlev)
@@ -923,7 +923,7 @@
         ! properties (noting that any gas scattering will have an
         ! asymmetry factor of zero)
-        do jlev = istartlev,iendlev
+        DO jlev = istartlev,iendlev
           if (od_sw_aerosol(1,jlev) > 0.0_jprb) then
-            do jg = 1,config%n_g_sw
+            DO jg = 1,config%n_g_sw
               iband = config%i_band_from_reordered_g_sw(jg)
               local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband,jlev)
@@ -961,9 +961,9 @@
 
         ! Loop over position
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
 ! Added for DWD (2020)
 !NEC$ forced_collapse
-          do jlev = istartlev,iendlev
-            do jb = 1,config%n_bands_lw
+          DO jlev = istartlev,iendlev
+            DO jb = 1,config%n_bands_lw
               od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol)
               scat_lw_aerosol(jb,jlev) = aerosol%ssa_lw(jb,jlev,jcol) * od_lw_aerosol(jb,jlev)
@@ -974,6 +974,6 @@
             end do
           end do
-          do jlev = istartlev,iendlev
-            do jg = 1,config%n_g_lw
+          DO jlev = istartlev,iendlev
+            DO jg = 1,config%n_g_lw
               iband = config%i_band_from_reordered_g_lw(jg)
               if (od_lw_aerosol(iband,jlev) > 0.0_jprb) then
@@ -995,18 +995,18 @@
 
         ! Loop over position
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
 ! Added for DWD (2020)
 !NEC$ forced_collapse
-          do jlev = istartlev,iendlev
+          DO jlev = istartlev,iendlev
             ! If aerosol longwave scattering is not included then we
             ! weight the optical depth by the single scattering
             ! co-albedo
-            do jb = 1, config%n_bands_lw
+            DO jb = 1, config%n_bands_lw
               od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol) &
                  &  * (1.0_jprb - aerosol%ssa_lw(jb,jlev,jcol))
             end do
           end do
-          do jlev = istartlev,iendlev
-            do jg = 1,config%n_g_lw
+          DO jlev = istartlev,iendlev
+            DO jg = 1,config%n_g_lw
               od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
                    &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev)
@@ -1120,5 +1120,5 @@
     if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_extinction',0,hook_handle)
 
-    do jtype = 1,config%n_aerosol_types
+    DO jtype = 1,config%n_aerosol_types
       if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
         write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
@@ -1138,10 +1138,10 @@
 
     ! Loop over position
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       ext = 0.0_jprb
       ! Get relative-humidity index
       irh = ao%calc_rh_index(relative_humidity(jcol))
       ! Add extinction coefficients from each aerosol type
-      do jtype = 1,config%n_aerosol_types
+      DO jtype = 1,config%n_aerosol_types
         if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
           ext = ext + mixing_ratio(jcol,jtype) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90	(revision 5158)
@@ -619,5 +619,5 @@
     iend   = ubound(itypes,1)
 
-    do jtype = istart, iend
+    DO jtype = istart, iend
       if (itypes(jtype) == 0) then
         call this%set_empty_type(jtype)
@@ -655,5 +655,5 @@
     else
       calc_rh_index = 1
-      do while (rh > this%rh_lower(calc_rh_index + 1))
+      DO while (rh > this%rh_lower(calc_rh_index + 1))
         calc_rh_index = calc_rh_index + 1
       end do
@@ -682,5 +682,5 @@
     end if
 
-    do jtype = 1,size(i_type_map)
+    DO jtype = 1,size(i_type_map)
       if (i_type_map(jtype) > 0) then
         write(nulout,'(i4,a,a)') jtype, ' -> hydrophobic type ', &
@@ -714,5 +714,5 @@
 
     ! Find index of first character
-    do while (i_line_current < iline)
+    DO while (i_line_current < iline)
       i_start_new = scan(str(istart:iend), new_line(' '))
       if (i_start_new == 0) then
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90	(revision 5158)
@@ -165,5 +165,5 @@
     
     ! Loop over hydrophilic types
-    do ja = 1,size(this%bin_philic)
+    DO ja = 1,size(this%bin_philic)
       ! Check if we have a match
       if (to_string(this%code_philic(:,ja)) == code_str &
@@ -176,5 +176,5 @@
     end do
     ! Repeat for the hydrophobic types
-    do ja = 1,size(this%bin_phobic)
+    DO ja = 1,size(this%bin_phobic)
       if (to_string(this%code_phobic(:,ja)) == code_str &
            &  .and. trim(to_string(this%optical_model_phobic(:,ja))) &
@@ -255,5 +255,5 @@
     if (lhydrophilic) then
       ! Loop over hydrophilic aerosol types
-      do ja = 1,size(this%bin_philic)
+      DO ja = 1,size(this%bin_philic)
         current_score = 0
         if (to_string(this%code_philic(:,ja)) == code_str) then
@@ -306,5 +306,5 @@
     else
       ! Loop over hydrophobic aerosol types
-      do ja = 1,size(this%bin_phobic)
+      DO ja = 1,size(this%bin_phobic)
         current_score = 0
         if (to_string(this%code_phobic(:,ja)) == code_str) then
@@ -376,5 +376,5 @@
     character(len=size(arr)) :: str
     integer :: jc
-    do jc = 1,size(arr)
+    DO jc = 1,size(arr)
       if (ichar(arr(jc)) == 0) then
         ! Replace NULL character with a space
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud.F90	(revision 5158)
@@ -251,5 +251,5 @@
       ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
       ! don't take the logarithm of the first pressure in each column.
-      do jcol = i1,i2
+      DO jcol = i1,i2
         this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length) &
              &                            * thermodynamics%temperature_hl(jcol,2) &
@@ -258,6 +258,6 @@
       end do
 
-      do jlev = 2,nlev-1
-        do jcol = i1,i2
+      DO jlev = 2,nlev-1
+        DO jcol = i1,i2
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -271,6 +271,6 @@
        ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
        ! don't take the logarithm of the last pressure in each column.
-      do jlev = 1,nlev-2
-        do jcol = i1,i2
+      DO jlev = 1,nlev-2
+        DO jcol = i1,i2
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -280,5 +280,5 @@
       end do
 
-      do jcol = i1,i2
+      DO jcol = i1,i2
         this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length) &
             &                            * thermodynamics%temperature_hl(jcol,nlev) &
@@ -340,5 +340,5 @@
       ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
       ! don't take the logarithm of the first pressure in each column.
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol)) &
              &                            * thermodynamics%temperature_hl(jcol,2) &
@@ -347,6 +347,6 @@
       end do
 
-      do jlev = 2,nlev-1
-        do jcol = istartcol,iendcol
+      DO jlev = 2,nlev-1
+        DO jcol = istartcol,iendcol
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -360,6 +360,6 @@
        ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
        ! don't take the logarithm of the last pressure in each column.
-      do jlev = 1,nlev-2
-        do jcol = istartcol,iendcol
+      DO jlev = 1,nlev-2
+        DO jcol = istartcol,iendcol
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -369,5 +369,5 @@
       end do
 
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol)) &
             &                            * thermodynamics%temperature_hl(jcol,nlev) &
@@ -423,5 +423,5 @@
       ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
       ! don't take the logarithm of the first pressure in each column.
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol,1)) &
              &                            * thermodynamics%temperature_hl(jcol,2) &
@@ -430,6 +430,6 @@
       end do
 
-      do jlev = 2,nlev-1
-        do jcol = istartcol,iendcol
+      DO jlev = 2,nlev-1
+        DO jcol = istartcol,iendcol
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -443,6 +443,6 @@
        ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
        ! don't take the logarithm of the last pressure in each column.
-      do jlev = 1,nlev-2
-        do jcol = istartcol,iendcol
+      DO jlev = 1,nlev-2
+        DO jcol = istartcol,iendcol
           this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) &
               &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
@@ -452,5 +452,5 @@
       end do
 
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
       ! AI ATTENTION a verifier decorrelation_length(jcol,nlev-1) ou nlev
         this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol,nlev-1)) &
@@ -663,5 +663,5 @@
     end if
 
-    do jcol = i1,i2
+    DO jcol = i1,i2
       eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
            &  * (0.5_jprb / pressure_hl(jcol,isurf))
@@ -761,5 +761,5 @@
     end if
 
-    do jcol = i1,i2
+    DO jcol = i1,i2
       eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
            &  * (0.5_jprb / pressure_hl(jcol,isurf))
@@ -804,14 +804,14 @@
     ntype = size(this%mixing_ratio,3)
     
-    do jlev = 1,nlev
-      do jcol = istartcol,iendcol
+    DO jlev = 1,nlev
+      DO jcol = istartcol,iendcol
         sum_mixing_ratio(jcol) = 0.0_jprb
       end do
-      do jh = 1, ntype
-        do jcol = istartcol,iendcol
+      DO jh = 1, ntype
+        DO jcol = istartcol,iendcol
           sum_mixing_ratio(jcol) = sum_mixing_ratio(jcol) + this%mixing_ratio(jcol,jlev,jh)
         end do
       end do
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         if (this%fraction(jcol,jlev)        < cloud_fraction_threshold &
              &  .or. sum_mixing_ratio(jcol) < cloud_mixing_ratio_threshold) then
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90	(revision 5158)
@@ -159,5 +159,5 @@
     cum_product = 1.0_jprb - frac(1)
     cum_cloud_cover(1) = frac(1)
-    do jlev = 1,nlev-1
+    DO jlev = 1,nlev-1
       ! Compute the combined cloud cover of layers jlev and jlev+1
       pair_cloud_cover(jlev) = max(frac(jlev),frac(jlev+1))
@@ -242,5 +242,5 @@
     cum_product = 1.0_jprb - frac(1)
     cum_cloud_cover(1) = frac(1)
-    do jlev = 1,nlev-1
+    DO jlev = 1,nlev-1
       ! Convert to "alpha" overlap parameter if necessary
       if (do_overlap_conversion) then
@@ -258,5 +258,5 @@
 #ifdef DWD_VECTOR_OPTIMIZATIONS
     end do
-    do jlev = 1,nlev-1
+    DO jlev = 1,nlev-1
 #endif
       if (frac(jlev) >= MaxCloudFrac) then
@@ -385,5 +385,5 @@
     jlev = 1
     nobj = 0
-    do while (jlev <= nlev)
+    DO while (jlev <= nlev)
       if (frac(jlev) > min_frac) then
         ! Starting a new object: store its top
@@ -392,5 +392,5 @@
         ! Find its maximum cloud fraction
         jlev = jlev + 1
-        do while (jlev <= nlev) 
+        DO while (jlev <= nlev)
           if (frac(jlev) < frac(jlev-1)) then
             exit
@@ -400,5 +400,5 @@
         i_max_obj(nobj) = jlev - 1
         ! Find its base
-        do while (jlev <= nlev)
+        DO while (jlev <= nlev)
           if (frac(jlev) > frac(jlev-1) .or. frac(jlev) <= min_frac) then
             exit
@@ -435,5 +435,5 @@
 
         ! Compute the combined cloud cover of pairs of layers
-        do jlev = 1,nlev-1
+        DO jlev = 1,nlev-1
           pair_cloud_cover(jlev) &
                &  = overlap_param(jlev)*max(frac(jlev),frac(jlev+1)) &
@@ -444,5 +444,5 @@
         ! adjacent objects as the product of the layerwise overlap
         ! parameters between their layers of maximum cloud fraction
-        do jobj = 1,nobj-1
+        DO jobj = 1,nobj-1
           alpha_obj(jobj) &
                &  = product(overlap_param(i_max_obj(jobj):i_max_obj(jobj+1)-1))
@@ -456,5 +456,5 @@
              &                     frac(1:nlev-1), frac(2:nlev))
         ! Compute the combined cloud cover of pairs of layers
-        do jlev = 1,nlev-1
+        DO jlev = 1,nlev-1
           pair_cloud_cover(jlev) &
                &  = overlap_alpha(jlev)*max(frac(jlev),frac(jlev+1)) &
@@ -465,5 +465,5 @@
         ! adjacent objects as the product of the layerwise overlap
         ! parameters between their layers of maximum cloud fraction
-        do jobj = 1,nobj-1
+        DO jobj = 1,nobj-1
           alpha_obj(jobj) &
                &  = product(overlap_alpha(i_max_obj(jobj):i_max_obj(jobj+1)-1))
@@ -476,7 +476,7 @@
       ! of each object: this will later be converted to the cumulative
       ! cloud cover working down from TOA
-      do jobj = 1,nobj
+      DO jobj = 1,nobj
         cum_cloud_cover(i_top_obj(jobj)) = frac(i_top_obj(jobj))
-        do jlev = i_top_obj(jobj), i_base_obj(jobj)-1
+        DO jlev = i_top_obj(jobj), i_base_obj(jobj)-1
           if (frac(jlev) >= MaxCloudFrac) then
             ! Cloud cover has reached one
@@ -496,5 +496,5 @@
       ! Sequentially combine objects until there is only one left
       ! covering the full vertical extent of clouds in the profile
-      do while (nobj > 1)
+      DO while (nobj > 1)
         ! Find the most correlated adjacent pair of objects
         alpha_max = 0.0_jprb
@@ -506,5 +506,5 @@
 
         jobj = 1
-        do while (jobj < nobj)
+        DO while (jobj < nobj)
           if (alpha_obj(jobj) > alpha_max) then
             alpha_max = alpha_obj(jobj)
@@ -533,5 +533,5 @@
         ! Scale the combined cloud cover of the lower object to
         ! account for its overlap with the upper object
-        do jlev = i_top_obj(iobj2),i_base_obj(iobj2)
+        DO jlev = i_top_obj(iobj2),i_base_obj(iobj2)
           cum_cloud_cover(jlev) = cum_cloud_cover(i_base_obj(iobj1)) &
                +  cum_cloud_cover(jlev) * scaling
@@ -554,5 +554,5 @@
       ! Ensure that the combined cloud cover of pairs of layers is
       ! consistent with the overhang
-      do jlev = 1,nlev-1
+      DO jlev = 1,nlev-1
         pair_cloud_cover(jlev) = max(pair_cloud_cover(jlev), &
              &     frac(jlev)+cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90	(revision 5158)
@@ -152,5 +152,5 @@
 
     total_cloud_cover = cum_cloud_cover(nlev);
-    do jlev = 1,nlev-1
+    DO jlev = 1,nlev-1
       overhang(jlev) = cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev)
     end do
@@ -166,10 +166,10 @@
       ! Find range of cloudy layers
       jlev = 1
-      do while (frac(jlev) <= 0.0_jprb) 
+      DO while (frac(jlev) <= 0.0_jprb)
         jlev = jlev + 1
       end do
       ibegin = jlev
       iend = jlev
-      do jlev = jlev+1,nlev
+      DO jlev = jlev+1,nlev
         if (frac(jlev) > 0.0_jprb) then
           iend = jlev
@@ -180,5 +180,5 @@
       overlap_param_inhom = overlap_param
 
-      do jlev = ibegin,iend-1
+      DO jlev = ibegin,iend-1
         if (overlap_param(jlev) > 0.0_jprb) then
           overlap_param_inhom(jlev) &
@@ -208,10 +208,10 @@
         
         ! Loop over ng columns
-        do jg = 1,ng
+        DO jg = 1,ng
           ! Find the cloud top height corresponding to the current
           ! random number, and store in itrigger
           trigger = rand_top(jg) * total_cloud_cover
           jlev = ibegin
-          do while (trigger > cum_cloud_cover(jlev) .and. jlev < iend)
+          DO while (trigger > cum_cloud_cover(jlev) .and. jlev < iend)
             jlev = jlev + 1
           end do
@@ -328,5 +328,5 @@
     ! Loop from the layer below the local cloud top down to the
     ! bottom-most cloudy layer
-    do jlev = itrigger+1,iend+1
+    DO jlev = itrigger+1,iend+1
       do_fill_od_scaling = .false.
       if (jlev <= iend) then
@@ -367,5 +367,5 @@
 
         ! Loop through the sequence of cloudy layers
-        do jcloud = 2,n_layers_to_scale
+        DO jcloud = 2,n_layers_to_scale
           ! Use second random number, and inhomogeneity overlap
           ! parameter, to decide whether the first random number
@@ -463,5 +463,5 @@
     ! Loop from the layer below the local cloud top down to the
     ! bottom-most cloudy layer
-    do jlev = itrigger+1,iend
+    DO jlev = itrigger+1,iend
       iy = iy+1
       if (is_cloudy(jlev-1)) then
@@ -494,5 +494,5 @@
         
     ! Loop through the sequence of cloudy layers
-    do jcloud = 2,n_layers_to_scale
+    DO jcloud = 2,n_layers_to_scale
       ! Use second random number, and inhomogeneity overlap
       ! parameter, to decide whether the first random number
@@ -679,5 +679,5 @@
 
     ! Loop down through layers starting at the first cloudy layer
-    do jlev = ibegin,iend
+    DO jlev = ibegin,iend
 
       if (is_any_cloud(jlev)) then
@@ -685,5 +685,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-        do jg = 1,ng
+        DO jg = 1,ng
           ! The intention is that all these operations are vectorizable,
           ! since all are vector operations on vectors of length ng...
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90	(revision 5158)
@@ -288,6 +288,6 @@
       end if
 
-      do jlev = 1,nlev
-        do jcol = istartcol,iendcol
+      DO jlev = 1,nlev
+        DO jcol = istartcol,iendcol
           ! Only do anything if cloud is present (assume that
           ! cloud%crop_cloud_fraction has already been called)
@@ -459,5 +459,5 @@
   ! Added for DWD (2020)
   !NEC$ shortloop
-              do jb = 1, config%n_bands_lw
+              DO jb = 1, config%n_bands_lw
                 od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) + od_lw_ice(jb)
                 if (scat_od_lw_liq(jb)+scat_od_lw_ice(jb) > 0.0_jprb) then
@@ -477,5 +477,5 @@
   ! Added for DWD (2020)
   !NEC$ shortloop
-              do jb = 1, config%n_bands_lw
+              DO jb = 1, config%n_bands_lw
                 od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) - scat_od_lw_liq(jb) &
                       &                   + od_lw_ice(jb) - scat_od_lw_ice(jb)
@@ -484,5 +484,5 @@
   ! Added for DWD (2020)
   !NEC$ shortloop
-            do jb = 1, config%n_bands_sw
+            DO jb = 1, config%n_bands_sw
               od_sw_cloud(jb,jlev,jcol) = od_sw_liq(jb) + od_sw_ice(jb)
               g_sw_cloud(jb,jlev,jcol) = (g_sw_liq(jb) * scat_od_sw_liq(jb) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_lw.F90	(revision 5158)
@@ -95,9 +95,9 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
 
       ! Compute the reflectance and transmittance of all layers,
       ! neglecting clouds
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (config%do_lw_aerosol_scattering) then
           ! Scattering case: first compute clear-sky reflectance,
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_sw.F90	(revision 5158)
@@ -103,5 +103,5 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       ! Only perform calculation if sun above the horizon
       if (single_level%cos_sza(jcol) > 0.0_jprb) then
@@ -114,5 +114,5 @@
           ! Delta-Eddington scaling has already been performed to the
           ! aerosol part of od, ssa and g
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             call calc_two_stream_gammas_sw(ng, cos_sza, &
                  &  ssa(:,jlev,jcol), g(:,jlev,jcol), &
@@ -128,5 +128,5 @@
         else
           ! Apply delta-Eddington scaling to the aerosol-gas mixture
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             od_total  =  od(:,jlev,jcol)
             ssa_total = ssa(:,jlev,jcol)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_config.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_config.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_config.F90	(revision 5158)
@@ -1657,5 +1657,5 @@
     ! contains the weights of interest.  We now find the non-zero weights
     nweights = 0
-    do jband = 1,size(mapping,2)
+    DO jband = 1,size(mapping,2)
       if (mapping(2,jband) > 0.0_jprb) then
         nweights = nweights+1
@@ -1674,9 +1674,9 @@
            &  wavenumber2, ' cm-1):'
       if (this%do_cloud_aerosol_per_sw_g_point) then
-        do jband = 1, nweights
+        DO jband = 1, nweights
           write(nulout, '(a,i0,a,f8.4)') '  Shortwave g point ', iband(jband), ': ', weight(jband)
         end do
       else
-        do jband = 1, nweights
+        DO jband = 1, nweights
           wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband))
           wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband))
@@ -1728,5 +1728,5 @@
     allocate(diag_ind(ninterval+2))
     diag_ind = 0
-    do jint = 1,ninterval+2
+    DO jint = 1,ninterval+2
       diag_ind(jint) = jint
     end do
@@ -1894,5 +1894,5 @@
     ! Count the number of albedo/emissivity intervals
     ninterval = 0
-    do jint = 1,NMaxAlbedoIntervals
+    DO jint = 1,NMaxAlbedoIntervals
       if (this%i_sw_albedo_index(jint) > 0) then
         ninterval = jint
@@ -1948,5 +1948,5 @@
         write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), &
              &  ' shortwave intervals to albedo intervals:'
-        do jband = 1,size(this%i_albedo_from_band_sw)
+        DO jband = 1,size(this%i_albedo_from_band_sw)
           write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband)
         end do
@@ -1972,5 +1972,5 @@
     ! Count the number of albedo/emissivity intervals
     ninterval = 0
-    do jint = 1,NMaxAlbedoIntervals
+    DO jint = 1,NMaxAlbedoIntervals
       if (this%i_lw_emiss_index(jint) > 0) then
         ninterval = jint
@@ -2026,5 +2026,5 @@
         write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), &
              &  ' longwave intervals to emissivity intervals:'
-        do jband = 1,size(this%i_emiss_from_band_lw)
+        DO jband = 1,size(this%i_emiss_from_band_lw)
           write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband)
         end do
@@ -2056,5 +2056,5 @@
       is_not_found = .true.
 
-      do jc = 0,size(enum_str)-1
+      DO jc = 0,size(enum_str)-1
         if (trim(str) == trim(enum_str(jc))) then
           icode = jc
@@ -2066,5 +2066,5 @@
         write(nulerr,'(a,a,a,a,a)',advance='no') '*** Error: ', trim(var_name), &
              &  ' must be one of: "', enum_str(0), '"'
-        do jc = 1,size(enum_str)-1
+        DO jc = 1,size(enum_str)-1
           write(nulerr,'(a,a,a)',advance='no') ', "', trim(enum_str(jc)), '"'
         end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_delta_eddington.h
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_delta_eddington.h	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_delta_eddington.h	(revision 5158)
@@ -86,5 +86,5 @@
   integer :: j
 
-  do j = 1,ng
+  DO j = 1,ng
     g            = scat_od_g(j) / max(scat_od(j), 1.0e-24)
     f            = g*g
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd.F90	(revision 5158)
@@ -209,5 +209,5 @@
     istart = 1
     this%i_gas_mapping = 0
-    do jgas = 1, this%ngas
+    DO jgas = 1, this%ngas
       if (jgas < this%ngas) then
         inext = istart + scan(constituent_id(istart:nchar), ' ')
@@ -217,5 +217,5 @@
       ! Find gas code
       i_gas_code = 0
-      do jjgas = 1, NMaxGases
+      DO jjgas = 1, NMaxGases
         if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
           i_gas_code = jjgas
@@ -264,5 +264,5 @@
          &  this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
     write(nulout, '(a)') '  Gases:'
-    do jgas = 1,this%ngas
+    DO jgas = 1,this%ngas
       if (this%single_gas(jgas)%i_gas_code > 0) then
         write(nulout, '(a,a,a)', advance='no') '    ', &
@@ -373,6 +373,6 @@
 
     ! Interpolate input SSI to regular wavenumber grid
-    do jwav_grid = 1,nwav_grid
-      do jwav = 1,nwav-1
+    DO jwav_grid = 1,nwav_grid
+      DO jwav = 1,nwav-1
         if (wavenumber(jwav) < wavenumber_grid(jwav_grid) &
              &  .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then
@@ -426,5 +426,5 @@
            &  ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance'
       iband = 0
-      do jg = 1,this%ng
+      DO jg = 1,this%ng
         if (this%spectral_def%i_band_number(jg) > iband) then
           iband = this%spectral_def%i_band_number(jg)
@@ -516,9 +516,9 @@
     global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
 
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
 
       log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
 
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         ! Find interpolation points in pressure
         pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
@@ -546,5 +546,5 @@
       optical_depth_fl(:,:,jcol) = 0.0_jprb
       
-      do jgas = 1,this%ngas
+      DO jgas = 1,this%ngas
 
         associate (single_gas => this%single_gas(jgas))
@@ -561,5 +561,5 @@
             end if
             
-            do jlev = 1,nlev
+            DO jlev = 1,nlev
               optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
                    &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
@@ -581,5 +581,5 @@
             end if
             
-            do jlev = 1,nlev
+            DO jlev = 1,nlev
               optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
                    &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
@@ -592,5 +592,5 @@
             ! Composite gases
             molar_abs => this%single_gas(jgas)%molar_abs
-            do jlev = 1,nlev
+            DO jlev = 1,nlev
               optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
                    &  + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
@@ -611,5 +611,5 @@
             molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
             mole_frac1 = exp(single_gas%log_mole_frac1)
-            do jlev = 1,nlev
+            DO jlev = 1,nlev
               ! Take care of mole_fraction == 0
               log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode)*scaling, mole_frac1))
@@ -628,5 +628,5 @@
               !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
               !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
-            do jlev = 1,nlev
+            DO jlev = 1,nlev
               optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
                    &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode) * scaling) * ( &
@@ -651,5 +651,5 @@
       ! Rayleigh scattering
       if (this%is_sw .and. present(rayleigh_od_fl)) then
-        do jlev = 1,nlev
+        DO jlev = 1,nlev
           rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
                &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
@@ -729,6 +729,6 @@
     od_fl(:,:,:) = 0.0_jprb
 
-    do jlev = 1,nlev
-      do jcol = istartcol,iendcol
+    DO jlev = 1,nlev
+      DO jcol = istartcol,iendcol
 
         log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1)))
@@ -758,5 +758,5 @@
     end do
 
-    do jgas = 1,this%ngas
+    DO jgas = 1,this%ngas
 
       associate (single_gas => this%single_gas(jgas))
@@ -768,7 +768,7 @@
           molar_abs => this%single_gas(jgas)%molar_abs
 
-          do jlev = 1,nlev
-            do jg = 1, this%ng
-              do jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jg = 1, this%ng
+              DO jcol = istartcol,iendcol
                 multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)
 
@@ -785,7 +785,7 @@
           molar_abs => this%single_gas(jgas)%molar_abs
 
-          do jlev = 1,nlev
-            do jg = 1, this%ng
-              do jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jg = 1, this%ng
+              DO jcol = istartcol,iendcol
                 multiplier = simple_multiplier(jcol,jlev)  * (mole_fraction_fl(jcol,jlev,igascode) &
                    &                            - single_gas%reference_mole_frac)
@@ -804,7 +804,7 @@
           molar_abs => this%single_gas(jgas)%molar_abs
 
-          do jlev = 1,nlev
-            do jg = 1, this%ng
-              do jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jg = 1, this%ng
+              DO jcol = istartcol,iendcol
                 od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
                    &        + (simple_multiplier(jcol,jlev)*tw1(jcol,jlev)) *   &
@@ -823,6 +823,6 @@
             mole_frac1 = exp(single_gas%log_mole_frac1)
 
-            do jlev = 1,nlev
-              do jcol = istartcol,iendcol
+            DO jlev = 1,nlev
+              DO jcol = istartcol,iendcol
                 ! Take care of mole_fraction == 0
                 log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
@@ -835,8 +835,8 @@
             end do
 
-          do jlev = 1,nlev
-            do jg = 1, this%ng
+          DO jlev = 1,nlev
+            DO jg = 1, this%ng
 !NEC$ select_vector
-              do jcol = istartcol,iendcol
+              DO jcol = istartcol,iendcol
 
                 od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
@@ -866,7 +866,7 @@
 
       ! Ensure the optical depth is not negative
-      do jcol = istartcol,iendcol
-        do jlev = 1,nlev
-          do jg = 1, this%ng
+      DO jcol = istartcol,iendcol
+        DO jlev = 1,nlev
+          DO jg = 1, this%ng
             optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev))
           end do
@@ -876,7 +876,7 @@
       ! Rayleigh scattering
       if (this%is_sw .and. present(rayleigh_od_fl)) then
-        do jcol = istartcol,iendcol
-          do jlev = 1,nlev
-            do jg = 1, this%ng
+        DO jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jg = 1, this%ng
               rayleigh_od_fl(jg,jlev,jcol) = global_multiplier &
                &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg)
@@ -908,5 +908,5 @@
     integer    :: it1, jt
 
-    do jt = 1,nt
+    DO jt = 1,nt
       tindex1 = (temperature(jt) - this%temperature1_planck) &
            &   * (1.0_jprb / this%d_temperature_planck)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90	(revision 5158)
@@ -272,7 +272,7 @@
       ! rayleigh optical depth: convert to total optical depth and
       ! single-scattering albedo
-      do jcol = istartcol,iendcol
-        do jlev = 1, nlev
-          do jg = 1, config%n_g_sw
+      DO jcol = istartcol,iendcol
+        DO jlev = 1, nlev
+          DO jg = 1, config%n_g_sw
             od_sw(jg,jlev,jcol)  = od_sw(jg,jlev,jcol) + ssa_sw(jg,jlev,jcol)
             ssa_sw(jg,jlev,jcol) = ssa_sw(jg,jlev,jcol) / od_sw(jg,jlev,jcol)
@@ -308,5 +308,5 @@
 
       ! Calculate the Planck function for each g point
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         call config%gas_optics_lw%calc_planck_function(nlev+1, &
              &  thermodynamics%temperature_hl(jcol,:), planck_hl(:,:,jcol))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_flux.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_flux.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_flux.F90	(revision 5158)
@@ -423,5 +423,5 @@
              &               config%i_band_from_reordered_g_sw, &
              &               this%sw_dn_surf_band, istartcol, iendcol)
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           this%sw_dn_surf_band(:,jcol) &
                &  = this%sw_dn_surf_band(:,jcol) &
@@ -429,5 +429,5 @@
         end do
       else
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
                &           config%i_band_from_reordered_g_sw, &
@@ -450,5 +450,5 @@
                &               config%i_band_from_reordered_g_sw, &
                &               this%sw_dn_surf_clear_band, istartcol, iendcol)
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             this%sw_dn_surf_clear_band(:,jcol) &
                  &  = this%sw_dn_surf_clear_band(:,jcol) &
@@ -456,5 +456,5 @@
           end do
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
                  &           config%i_band_from_reordered_g_sw, &
@@ -486,5 +486,5 @@
                &               this%sw_dn_diffuse_surf_canopy, istartcol, iendcol)
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
                  &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
@@ -502,6 +502,6 @@
         this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
         this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb
-        do jband = 1,config%n_bands_sw
-          do jalbedoband = 1,nalbedoband
+        DO jband = 1,config%n_bands_sw
+          DO jalbedoband = 1,nalbedoband
             if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
               ! Initially, "diffuse" is actually "total"
@@ -534,5 +534,5 @@
                &               this%lw_dn_surf_canopy, istartcol, iendcol)
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%lw_dn_surf_g(:,jcol), &
                  &           config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
@@ -548,5 +548,5 @@
                &               lw_dn_surf_band, istartcol, iendcol)
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%lw_dn_surf_g(:,jcol), &
                  &           config%i_band_from_reordered_g_lw, &
@@ -556,6 +556,6 @@
         nalbedoband = size(config%lw_emiss_weights,1)
         this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
-        do jband = 1,config%n_bands_lw
-          do jalbedoband = 1,nalbedoband
+        DO jband = 1,config%n_bands_lw
+          DO jalbedoband = 1,nalbedoband
             if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
               this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
@@ -602,5 +602,5 @@
              &               this%sw_up_toa_band, istartcol, iendcol)
       else
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           call indexed_sum(this%sw_dn_toa_g(:,jcol), &
                &           config%i_band_from_reordered_g_sw, &
@@ -618,5 +618,5 @@
                &               this%sw_up_toa_clear_band, istartcol, iendcol)
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%sw_up_toa_clear_g(:,jcol), &
                  &               config%i_band_from_reordered_g_sw, &
@@ -634,5 +634,5 @@
              &               this%lw_up_toa_band, istartcol, iendcol)
       else
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           call indexed_sum(this%lw_up_toa_g(:,jcol), &
                &           config%i_band_from_reordered_g_lw, &
@@ -647,5 +647,5 @@
                &               this%lw_up_toa_clear_band, istartcol, iendcol)
         else
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             call indexed_sum(this%lw_up_toa_clear_g(:,jcol), &
                  &               config%i_band_from_reordered_g_lw, &
@@ -753,5 +753,5 @@
     iend   = ubound(source,1)
 
-    do jg = istart, iend
+    DO jg = istart, iend
       ig = ind(jg)
       dest(ig) = dest(ig) + source(jg)
@@ -777,5 +777,5 @@
     iend   = ubound(source,1)
 
-    do jg = istart, iend
+    DO jg = istart, iend
       ig = ind(jg)
       dest(ig) = dest(ig) + source(jg)
@@ -797,7 +797,7 @@
     dest = 0.0
 
-    do jg = lbound(source,1), ubound(source,1)
+    DO jg = lbound(source,1), ubound(source,1)
       ig = ind(jg)
-      do jc = ist, iend
+      DO jc = ist, iend
         dest(ig,jc) = dest(ig,jc) + source(jg,jc)
       end do
@@ -820,6 +820,6 @@
     nlev   = size(source,2)
 
-    do jlev = 1,nlev
-      do jg = istart, iend
+    DO jlev = 1,nlev
+      DO jg = istart, iend
         ig = ind(jg)
         dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
@@ -846,6 +846,6 @@
     nlev   = size(source,2)
 
-    do jlev = 1,nlev
-      do jg = istart, iend
+    DO jlev = 1,nlev
+      DO jg = istart, iend
         ig = ind(jg)
         dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_gas.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_gas.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_gas.F90	(revision 5158)
@@ -209,6 +209,6 @@
     this%is_well_mixed(igas) = .false.
 
-    do jk = 1,this%nlev
-      do jc = i1,i2
+    DO jk = 1,this%nlev
+      DO jc = i1,i2
         this%mixing_ratio(jc,jk,igas) = mixing_ratio(jc-i1+1,jk)
       end do
@@ -294,6 +294,6 @@
     this%is_well_mixed(igas)           = .true.
 
-    do jk = 1,this%nlev
-      do jc = i1,i2
+    DO jk = 1,this%nlev
+      DO jc = i1,i2
         this%mixing_ratio(jc,jk,igas) = mixing_ratio
       end do
@@ -397,5 +397,5 @@
       end if
     else
-      do jg = 1,this%ntype
+      DO jg = 1,this%ntype
         call this%set_units(iunits, igas=this%icode(jg), scale_factor=new_sf)
       end do
@@ -416,5 +416,5 @@
     
     scaling = this%scale_factor
-    do jg = 1,NMaxGases
+    DO jg = 1,NMaxGases
       if (iunits == IMassMixingRatio .and. this%iunits(jg) == IVolumeMixingRatio) then
         scaling(jg) = scaling(jg) * GasMolarMass(jg) / AirMolarMass
@@ -481,5 +481,5 @@
       end if
     else
-      do jg = 1,this%ntype
+      DO jg = 1,this%ntype
         call this%assert_units(iunits, igas=this%icode(jg), scale_factor=sf, istatus=istatus)
       end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90	(revision 5158)
@@ -53,5 +53,5 @@
     ! Count number of cloud types
     config%n_cloud_types = 0
-    do jtype = 1,NMaxCloudTypes
+    DO jtype = 1,NMaxCloudTypes
       if (len_trim(config%cloud_type_name(jtype)) > 0) then
         config%n_cloud_types = jtype
@@ -83,5 +83,5 @@
 
     ! Load cloud optics data
-    do jtype = 1,config%n_cloud_types
+    DO jtype = 1,config%n_cloud_types
       if (config%cloud_type_name(jtype)(1:1) == '/') then
         file_name = trim(config%cloud_type_name(jtype))
@@ -191,5 +191,5 @@
 
     ! Loop over cloud types
-    do jtype = 1,config%n_cloud_types
+    DO jtype = 1,config%n_cloud_types
       ! Compute in-cloud water path
       if (config%is_homogeneous) then
@@ -238,6 +238,6 @@
     ! Scale the combined longwave optical properties
     if (config%do_lw_cloud_scattering) then
-      do jcol = istartcol, iendcol
-        do jlev = 1,nlev
+      DO jcol = istartcol, iendcol
+        DO jlev = 1,nlev
           if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
             ! Note that original cloud optics does not do
@@ -259,6 +259,6 @@
     if (config%do_sw) then
       if (.not. config%do_sw_delta_scaling_with_gases) then
-        do jcol = istartcol, iendcol
-          do jlev = 1,nlev
+        DO jcol = istartcol, iendcol
+          DO jlev = 1,nlev
             if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
               call delta_eddington_extensive(od_sw_cloud(:,jlev,jcol), &
@@ -269,9 +269,9 @@
       end if
 
-      do jcol = istartcol, iendcol
-        do jlev = 1,nlev
+      DO jcol = istartcol, iendcol
+        DO jlev = 1,nlev
           if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
             ! Scale to get asymmetry factor and single scattering albedo
-            do jg = 1, config%n_bands_sw
+            DO jg = 1, config%n_bands_sw
               g_sw_cloud(jg,jlev,jcol) = g_sw_cloud(jg,jlev,jcol) &
                  &  / max(ssa_sw_cloud(jg,jlev,jcol), 1.0e-15_jprb)
@@ -308,5 +308,5 @@
     if (lhook) call dr_hook('radiation_general_cloud_optics:save',0,hook_handle)
 
-    do jtype = 1,config%n_cloud_types
+    DO jtype = 1,config%n_cloud_types
       if (config%do_sw) then
         associate(co_sw => config%cloud_optics_sw(jtype))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90	(revision 5158)
@@ -285,6 +285,6 @@
 
     if (present(scat_od)) then
-      do jcol = 1,ncol
-        do jlev = 1,nlev
+      DO jcol = 1,ncol
+        DO jlev = 1,nlev
           if (cloud_fraction(jcol, jlev) > 0.0_jprb) then
             re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) &
@@ -293,5 +293,5 @@
             weight2 = re_index - ire
             weight1 = 1.0_jprb - weight2
-            do jg = 1, ng
+            DO jg = 1, ng
               od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(jg,ire) &
                  &                                +weight2*this%mass_ext(jg,ire+1))
@@ -310,6 +310,6 @@
     else
       ! No scattering: return the absorption optical depth
-      do jcol = 1,ncol
-        do jlev = 1,nlev
+      DO jcol = 1,ncol
+        DO jlev = 1,nlev
           if (water_path(jcol, jlev) > 0.0_jprb) then
             re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) &
@@ -404,5 +404,5 @@
 
     ! Define effective radius
-    do ire = 1,this%n_effective_radius
+    DO ire = 1,this%n_effective_radius
       effective_radius(ire) = this%effective_radius_0 + this%d_effective_radius*(ire-1)
     end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90	(revision 5158)
@@ -119,9 +119,9 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
 
       ! Is there any cloud in the profile?
       is_cloudy_profile = .false.
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
           is_cloudy_profile = .true.
@@ -135,5 +135,5 @@
       ! the clear-sky layers since these will be needed when we come
       ! to do the total-sky fluxes.
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (config%do_clear .or. cloud%fraction(jcol,jlev) &
              &                 < config%cloud_fraction_threshold) then
@@ -201,5 +201,5 @@
       ! now.
       if (is_cloudy_profile .or. .not. config%do_clear) then
-        do jlev = 1,nlev
+        DO jlev = 1,nlev
           ! Compute combined gas+aerosol+cloud optical properties;
           ! note that for clear layers, the reflectance and
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90	(revision 5158)
@@ -123,5 +123,5 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       ! Only perform calculation if sun above the horizon
       if (single_level%cos_sza(jcol) > 0.0_jprb) then
@@ -131,5 +131,5 @@
         ! Is there any cloud in the profile?
         is_cloudy_profile = .false.
-        do jlev = 1,nlev
+        DO jlev = 1,nlev
           if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
             is_cloudy_profile = .true.
@@ -146,5 +146,5 @@
           ! Delta-Eddington scaling has already been performed to the
           ! aerosol part of od, ssa and g
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             if (config%do_clear .or. cloud%fraction(jcol,jlev) &
                  &                 < config%cloud_fraction_threshold) then
@@ -164,5 +164,5 @@
         else
           ! Apply delta-Eddington scaling to the aerosol-gas mixture
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             if (config%do_clear .or. cloud%fraction(jcol,jlev) &
                  &                 < config%cloud_fraction_threshold) then
@@ -229,5 +229,5 @@
         ! fluxes now.
         if (is_cloudy_profile .or. .not. config%do_clear) then
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             ! Compute combined gas+aerosol+cloud optical properties;
             ! note that for clear layers, the reflectance and
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90	(revision 5158)
@@ -73,5 +73,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jb = 1, nb
+    DO jb = 1, nb
       od(jb) = iwp_gm_2 * (coeff(jb,1) + coeff(jb,2) * inv_de_um)
       scat_od(jb) = od(jb) * (1.0_jprb - (coeff(jb,3) + de_um*(coeff(jb,4) &
@@ -124,5 +124,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jb = 1, nb
+    DO jb = 1, nb
       od(jb) = iwp_gm_2 * (coeff(jb,1) + inv_de_um*(coeff(jb,2) &
          &  + inv_de_um*coeff(jb,3)))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90	(revision 5158)
@@ -386,6 +386,6 @@
     ZONEMINUS_ARRAY = ZONEMINUS
 
-    do jlev=1,nlev
-      do jcol= istartcol,iendcol
+    DO jlev=1,nlev
+      DO jcol= istartcol,iendcol
         pressure_fl(jcol,jlev) &
             &  = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) &
@@ -486,11 +486,11 @@
         ! arrays have dimensions in a different order to the inputs,
         ! so there is some inefficiency here.
-        do jgreorder = 1,config%n_g_lw
+        DO jgreorder = 1,config%n_g_lw
           iband = config%i_band_from_reordered_g_lw(jgreorder)
           ig = config%i_g_from_reordered_g_lw(jgreorder)
           
           ! Top-of-atmosphere half level
-          do jlev = 1,nlev
-            do jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jcol = istartcol,iendcol
               ! Some g points can return negative optical depths;
               ! specifically original g points 54-56 which causes
@@ -504,6 +504,6 @@
       else
         ! G points have not been reordered 
-        do jcol = istartcol,iendcol
-          do jlev = 1,nlev
+        DO jcol = istartcol,iendcol
+          DO jlev = 1,nlev
             ! Check for negative optical depth
             od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
@@ -544,6 +544,6 @@
       ! Scale the incoming solar per band, if requested
       if (config%use_spectral_solar_scaling) then
-        do jg = 1,JPGPT_SW
-          do jcol = istartcol,iendcol 
+        DO jg = 1,JPGPT_SW
+          DO jcol = istartcol,iendcol
             ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
                  &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
@@ -556,5 +556,5 @@
       ! ZINCSOL will be zero.
       if (present(incoming_sw)) then
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           if (single_level%cos_sza(jcol) > 0.0_jprb) then
 ! Added for DWD (2020)
@@ -570,8 +570,8 @@
 !      if (.true.) then
         ! Account for reordered g points
-        do jgreorder = 1,config%n_g_sw
+        DO jgreorder = 1,config%n_g_sw
           ig = config%i_g_from_reordered_g_sw(jgreorder)
-          do jlev = 1,nlev
-            do jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jcol = istartcol,iendcol
               ! Check for negative optical depth
               od_sw (jgreorder,nlev+1-jlev,jcol) &
@@ -587,7 +587,7 @@
       else
         ! G points have not been reordered
-        do jcol = istartcol,iendcol
-          do jlev = 1,nlev
-            do jg = 1,config%n_g_sw
+        DO jcol = istartcol,iendcol
+          DO jlev = 1,nlev
+            DO jg = 1,config%n_g_sw
               ! Check for negative optical depth
               od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
@@ -596,5 +596,5 @@
           end do
           if (present(incoming_sw)) then
-            do jg = 1,config%n_g_sw
+            DO jg = 1,config%n_g_sw
               incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
             end do
@@ -667,6 +667,6 @@
     ! lowest interpolation bound, and the fraction into interpolation
     ! interval
-    do jlev = 1,nlev+1
-      do jcol = istartcol,iendcol
+    DO jlev = 1,nlev+1
+      DO jcol = istartcol,iendcol
         temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset)
         if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then
@@ -687,7 +687,7 @@
 
       ! Calculate Planck functions per band
-      do jband = 1,config%n_bands_lw
+      DO jband = 1,config%n_bands_lw
         factor = zfluxfac * delwave(jband)
-        do jcol = istartcol,iendcol
+        DO jcol = istartcol,iendcol
           planck_store(jcol,jband) = factor &
                &  * (totplnk(ind(jcol),jband) &
@@ -706,5 +706,5 @@
           ! Top-of-atmosphere half level - note that PFRAC is on model
           ! levels not half levels
-          do jgreorder = 1,config%n_g_lw
+          DO jgreorder = 1,config%n_g_lw
             iband = config%i_band_from_reordered_g_lw(jgreorder)
             ig = config%i_g_from_reordered_g_lw(jgreorder)
@@ -713,5 +713,5 @@
           end do
         else
-          do jgreorder = 1,config%n_g_lw
+          DO jgreorder = 1,config%n_g_lw
             iband = config%i_band_from_reordered_g_lw(jgreorder)
             ig = config%i_g_from_reordered_g_lw(jgreorder)
@@ -726,14 +726,14 @@
           ! Top-of-atmosphere half level - note that PFRAC is on model
           ! levels not half levels
-          do jg = 1,config%n_g_lw
+          DO jg = 1,config%n_g_lw
             iband = config%i_band_from_g_lw(jg)
             planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev)
           end do
         else
-          do jg = 1,config%n_g_lw
+          DO jg = 1,config%n_g_lw
             iband = config%i_band_from_g_lw(jg)
             planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
           end do
-          do jcol = istartcol,iendcol
+          DO jcol = istartcol,iendcol
             planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
           end do
@@ -794,5 +794,5 @@
 
     ! Work out surface interpolations
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       Tsurf = temperature(jcol)
       if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then
@@ -813,7 +813,7 @@
 
     ! Calculate Planck functions per band
-    do jband = 1,config%n_bands_lw
+    DO jband = 1,config%n_bands_lw
       factor = zfluxfac * delwave(jband)
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         planck_store(jcol,jband) = factor &
              &  * (totplnk(ind(jcol),jband) &
@@ -829,5 +829,5 @@
       ! in pressure since the the functions above treat pressure
       ! decreasing with increasing index.
-      do jgreorder = 1,config%n_g_lw
+      DO jgreorder = 1,config%n_g_lw
         iband = config%i_band_from_reordered_g_lw(jgreorder)
         ig = config%i_g_from_reordered_g_lw(jgreorder)
@@ -836,5 +836,5 @@
     else
       ! G points have not been reordered 
-      do jg = 1,config%n_g_lw
+      DO jg = 1,config%n_g_lw
         iband = config%i_band_from_g_lw(jg)
         planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_liquid_optics_socrates.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_liquid_optics_socrates.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_liquid_optics_socrates.F90	(revision 5158)
@@ -65,5 +65,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jb = 1, nb
+    DO jb = 1, nb
       od(jb) = lwp * (coeff(jb,1) + re*(coeff(jb,2) + re*coeff(jb,3))) &
          &  / (1.0_jprb + re*(coeff(jb,4) + re*(coeff(jb,5) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_lw_derivatives.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_lw_derivatives.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_lw_derivatives.F90	(revision 5158)
@@ -74,5 +74,5 @@
     ! Move up through the atmosphere computing the derivatives at each
     ! half-level
-    do jlev = nlev,1,-1
+    DO jlev = nlev,1,-1
       lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev)
       lw_derivatives(icol,jlev) = sum(lw_derivatives_g)
@@ -122,5 +122,5 @@
     ! Move up through the atmosphere computing the derivatives at each
     ! half-level
-    do jlev = nlev,1,-1
+    DO jlev = nlev,1,-1
       lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev)
       lw_derivatives(icol,jlev) = (1.0_jprb - weight) * lw_derivatives(icol,jlev) &
@@ -177,5 +177,5 @@
     ! Move up through the atmosphere computing the derivatives at each
     ! half-level
-    do jlev = nlev,1,-1
+    DO jlev = nlev,1,-1
       ! Compute effect of overlap at half-level jlev+1, yielding
       ! derivatives just above that half-level
@@ -243,5 +243,5 @@
       ! multiplication inline
       
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
         ! Compute effect of overlap at half-level jlev+1, yielding
         ! derivatives just above that half-level
@@ -249,5 +249,5 @@
         
         associate(A=>u_matrix(:,:,jlev+1), b=>lw_deriv_below)
-          do jg = 1,ng   
+          DO jg = 1,ng
             ! Both inner and outer loop of the matrix loops j1 and j2 unrolled
             ! inner loop:        j2=1             j2=2             j2=3 
@@ -273,5 +273,5 @@
       ! Move up through the atmosphere computing the derivatives at each
       ! half-level
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
         ! Compute effect of overlap at half-level jlev+1, yielding
         ! derivatives just above that half-level
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_matrix.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_matrix.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_matrix.F90	(revision 5158)
@@ -91,6 +91,6 @@
       mat_x_vec(1:iend,1) = A(1:iend,1,1)*b(1:iend,1)
     else
-      do j1 = 1,m
-        do j2 = 1,m
+      DO j1 = 1,m
+        DO j2 = 1,m
           mat_x_vec(1:iend,j1) = mat_x_vec(1:iend,j1) &
                &               + A(1:iend,j1,j2)*b(1:iend,j2)
@@ -125,6 +125,6 @@
     singlemat_x_vec = 0.0_jprb
 
-    do j1 = 1,m
-      do j2 = 1,m
+    DO j1 = 1,m
+      DO j2 = 1,m
         singlemat_x_vec(1:iend,j1) = singlemat_x_vec(1:iend,j1) &
              &                    + A(j1,j2)*b(1:iend,j2)
@@ -176,7 +176,7 @@
       m2block = 2*mblock 
       ! Do the top-left (C, D, F, G)
-      do j2 = 1,m2block
-        do j1 = 1,m2block
-          do j3 = 1,m2block
+      DO j2 = 1,m2block
+        DO j1 = 1,m2block
+          DO j3 = 1,m2block
             mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
                  &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
@@ -184,8 +184,8 @@
         end do
       end do
-      do j2 = m2block+1,m
+      DO j2 = m2block+1,m
         ! Do the top-right (E & H)
-        do j1 = 1,m2block
-          do j3 = 1,m
+        DO j1 = 1,m2block
+          DO j3 = 1,m
             mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
                  &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
@@ -193,6 +193,6 @@
         end do
         ! Do the bottom-right (I)
-        do j1 = m2block+1,m
-          do j3 = m2block+1,m
+        DO j1 = m2block+1,m
+          DO j3 = m2block+1,m
             mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
                  &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
@@ -202,7 +202,7 @@
     else
       ! Ordinary dense matrix
-      do j2 = 1,m
-        do j1 = 1,m
-          do j3 = 1,m
+      DO j2 = 1,m
+        DO j1 = 1,m
+          DO j3 = 1,m
             mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
                  &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
@@ -238,7 +238,7 @@
     singlemat_x_mat = 0.0_jprb
 
-    do j2 = 1,m
-      do j1 = 1,m
-        do j3 = 1,m
+    DO j2 = 1,m
+      DO j1 = 1,m
+        DO j3 = 1,m
           singlemat_x_mat(1:iend,j1,j2) = singlemat_x_mat(1:iend,j1,j2) &
                &                        + A(j1,j3)*B(1:iend,j3,j2)
@@ -273,7 +273,7 @@
     mat_x_singlemat = 0.0_jprb
 
-    do j2 = 1,m
-      do j1 = 1,m
-        do j3 = 1,m
+    DO j2 = 1,m
+      DO j1 = 1,m
+        DO j3 = 1,m
           mat_x_singlemat(1:iend,j1,j2) = mat_x_singlemat(1:iend,j1,j2) &
                &                        + A(1:iend,j1,j3)*B(j3,j2)
@@ -311,5 +311,5 @@
 
     identity_minus_mat_x_mat = - identity_minus_mat_x_mat
-    do j = 1,m
+    DO j = 1,m
       identity_minus_mat_x_mat(1:iend,j,j) &
            &     = 1.0_jprb + identity_minus_mat_x_mat(1:iend,j,j)
@@ -337,6 +337,6 @@
     
     sparse_x_dense = 0.0_jprb
-    do j2 = 1,n2
-      do j1 = 1,n1
+    DO j2 = 1,n2
+      DO j1 = 1,n1
         if (sparse(j1,j2) /= 0.0_jprb) then
           sparse_x_dense(j1,:) = sparse_x_dense(j1,:) + sparse(j1,j2)*dense(j2,:)
@@ -376,10 +376,10 @@
       mblock = m/3
       m2block = 2*mblock
-      do j4 = 1,nrepeat
+      DO j4 = 1,nrepeat
         repeated_square = 0.0_jprb
         ! Do the top-left (C, D, F & G)
-        do j2 = 1,m2block
-          do j1 = 1,m2block
-            do j3 = 1,m2block
+        DO j2 = 1,m2block
+          DO j1 = 1,m2block
+            DO j3 = 1,m2block
               repeated_square(j1,j2) = repeated_square(j1,j2) &
                    &                 + A(j1,j3)*A(j3,j2)
@@ -387,8 +387,8 @@
           end do
         end do
-        do j2 = m2block+1, m
+        DO j2 = m2block+1, m
           ! Do the top-right (E & H)
-          do j1 = 1,m2block
-            do j3 = 1,m
+          DO j1 = 1,m2block
+            DO j3 = 1,m
               repeated_square(j1,j2) = repeated_square(j1,j2) &
                    &                 + A(j1,j3)*A(j3,j2)
@@ -396,6 +396,6 @@
           end do
           ! Do the bottom-right (I)
-          do j1 = m2block+1, m
-            do j3 = m2block+1,m
+          DO j1 = m2block+1, m
+            DO j3 = m2block+1,m
               repeated_square(j1,j2) = repeated_square(j1,j2) &
                    &                 + A(j1,j3)*A(j3,j2)
@@ -409,9 +409,9 @@
     else
       ! Ordinary dense matrix
-      do j4 = 1,nrepeat
+      DO j4 = 1,nrepeat
         repeated_square = 0.0_jprb
-        do j2 = 1,m
-          do j1 = 1,m
-            do j3 = 1,m
+        DO j2 = 1,m
+          DO j1 = 1,m
+            DO j3 = 1,m
               repeated_square(j1,j2) = repeated_square(j1,j2) &
                    &                 + A(j1,j3)*A(j3,j2)
@@ -549,5 +549,5 @@
     U33 = A(1:iend,3,3) - L31*A(1:iend,1,3) - L32*U23
 
-    do j = 1,3
+    DO j = 1,3
       ! Solve Ly = B(:,:,j) by forward substitution
       ! y1 = B(:,1,j)
@@ -649,15 +649,15 @@
     LU(1:iend,1:m,1:m) = A(1:iend,1:m,1:m)
 
-    do j2 = 1, m
-      do j1 = 1, j2-1
+    DO j2 = 1, m
+      DO j1 = 1, j2-1
         s = LU(1:iend,j1,j2)
-        do j3 = 1, j1-1
+        DO j3 = 1, j1-1
           s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2)
         end do
         LU(1:iend,j1,j2) = s
       end do
-      do j1 = j2, m
+      DO j1 = j2, m
         s = LU(1:iend,j1,j2)
-        do j3 = 1, j2-1
+        DO j3 = 1, j2-1
           s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2)
         end do
@@ -666,5 +666,5 @@
       if (j2 /= m) then
         s = 1.0_jprb / LU(1:iend,j2,j2)
-        do j1 = j2+1, m
+        DO j1 = j2+1, m
           LU(1:iend,j1,j2) = LU(1:iend,j1,j2) * s
         end do
@@ -691,12 +691,12 @@
 
     ! First solve Ly=b
-    do j2 = 2, m
-      do j1 = 1, j2-1
+    DO j2 = 2, m
+      DO j1 = 1, j2-1
         x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1)
       end do
     end do
     ! Now solve Ux=y
-    do j2 = m, 1, -1
-      do j1 = j2+1, m
+    DO j2 = m, 1, -1
+      DO j1 = j2+1, m
         x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1)
       end do
@@ -722,5 +722,5 @@
     call lu_factorization(n,iend,m,A,LU)
 
-    do j = 1, m
+    DO j = 1, m
       call lu_substitution(n,iend,m,LU,B(1:,1:m,j),X(1:iend,1:m,j))
 !      call lu_substitution(n,iend,m,LU,B(1:n,1:m,j),X(1:iend,1:m,j))
@@ -835,12 +835,12 @@
 
     ! Compute the 1-norms of A
-    do j3 = 1,m
+    DO j3 = 1,m
       sum_column(:) = 0.0_jprb
-      do j2 = 1,m
-        do j1 = 1,iend
+      DO j2 = 1,m
+        DO j1 = 1,iend
           sum_column(j1) = sum_column(j1) + abs(A(j1,j2,j3))
         end do
       end do
-      do j1 = 1,iend
+      DO j1 = 1,iend
         if (sum_column(j1) > normA(j1)) then
           normA(j1) = sum_column(j1)
@@ -861,6 +861,6 @@
     ! Scale the input matrices by a power of 2
     scaling = 2.0_jprb**(-expo)
-    do j3 = 1,m
-      do j2 = 1,m
+    DO j3 = 1,m
+      DO j2 = 1,m
         A(1:iend,j2,j3) = A(1:iend,j2,j3) * scaling
       end do
@@ -872,5 +872,5 @@
 
     V = c(8)*A6 + c(6)*A4 + c(4)*A2
-    do j3 = 1,m
+    DO j3 = 1,m
       V(:,j3,j3) = V(:,j3,j3) + c(2)
     end do
@@ -878,5 +878,5 @@
     V = c(7)*A6 + c(5)*A4 + c(3)*A2
     ! Add a multiple of the identity matrix
-    do j3 = 1,m
+    DO j3 = 1,m
       V(:,j3,j3) = V(:,j3,j3) + c(1)
     end do
@@ -887,10 +887,10 @@
 
     ! Add the identity matrix
-    do j3 = 1,m
+    DO j3 = 1,m
       A(1:iend,j3,j3) = A(1:iend,j3,j3) + 1.0_jprb
     end do
 
     ! Loop through the matrices
-    do j1 = 1,iend
+    DO j1 = 1,iend
       if (expo(j1) > 0) then
         ! Square matrix j1 expo(j1) times          
@@ -1016,6 +1016,6 @@
 
     ! Compute V * diag_rdivide_V
-    do j1 = 1,3
-      do j2 = 1,3
+    DO j1 = 1,3
+      DO j2 = 1,3
         R(1:iend,j2,j1) = V(1:iend,j2,1)*diag_rdivide_V(1:iend,1,j1) &
              &          + V(1:iend,j2,2)*diag_rdivide_V(1:iend,2,j1) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90	(revision 5158)
@@ -146,5 +146,5 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
 
       ! Clear-sky calculation
@@ -203,5 +203,5 @@
         is_clear_sky_layer = .true.
         i_cloud_top = nlev+1
-        do jlev = 1,nlev
+        DO jlev = 1,nlev
           ! Compute combined gas+aerosol+cloud optical properties
           if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
@@ -212,5 +212,5 @@
             end if
 
-            do jg = 1,ng
+            DO jg = 1,ng
               od_cloud_new(jg) = od_scaling(jg,jlev) &
                  &  * od_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol)
@@ -228,5 +228,5 @@
                 ! case that od_total > 0.0 and ssa_total > 0.0 but
                 ! od_total*ssa_total == 0 due to underflow
-                do jg = 1,ng
+                DO jg = 1,ng
                   if (od_total(jg) > 0.0_jprb) then
                     scat_od_total(jg) = ssa(jg,jlev,jcol)*od(jg,jlev,jcol) &
@@ -247,5 +247,5 @@
               else
 
-                do jg = 1,ng
+                DO jg = 1,ng
                   if (od_total(jg) > 0.0_jprb) then
                     scat_od = ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90	(revision 5158)
@@ -138,5 +138,5 @@
 
     ! Loop through columns
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       ! Only perform calculation if sun above the horizon
       if (single_level%cos_sza(jcol) > 0.0_jprb) then
@@ -155,5 +155,5 @@
         else
           ! Apply delta-Eddington scaling to the aerosol-gas mixture
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             od_total  =  od(:,jlev,jcol)
             ssa_total = ssa(:,jlev,jcol)
@@ -205,8 +205,8 @@
         if (total_cloud_cover >= config%cloud_fraction_threshold) then
           ! Total-sky calculation
-          do jlev = 1,nlev
+          DO jlev = 1,nlev
             ! Compute combined gas+aerosol+cloud optical properties
             if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
-              do jg = 1,ng
+              DO jg = 1,ng
                 od_cloud_new(jg) = od_scaling(jg,jlev) &
                    &  * od_cloud(config%i_band_from_reordered_g_sw(jg),jlev,jcol)
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_monochromatic.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_monochromatic.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_monochromatic.F90	(revision 5158)
@@ -164,5 +164,5 @@
     integer :: jlev
 
-    do jlev = 1,nlev
+    DO jlev = 1,nlev
       ! The fraction of the total optical depth in the current layer
       ! is proportional to the fraction of the mass of the atmosphere
@@ -192,5 +192,5 @@
              &  = StefanBoltzmann * single_level%skin_temperature(istartcol:iendcol)**4 &
              &  * single_level%lw_emissivity(istartcol:iendcol,1)
-        do jlev = 1,nlev+1
+        DO jlev = 1,nlev+1
           planck_hl(1,jlev,istartcol:iendcol) = StefanBoltzmann * thermodynamics%temperature_hl(istartcol:iendcol,jlev)**4
         end do
@@ -200,5 +200,5 @@
              &             single_level%skin_temperature(istartcol:iendcol)) &
              &  * single_level%lw_emissivity(istartcol:iendcol,1)
-        do jlev = 1,nlev+1
+        DO jlev = 1,nlev+1
           planck_hl(1,jlev,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, &
                &             thermodynamics%temperature_hl(istartcol:iendcol,jlev))
@@ -261,6 +261,6 @@
     ! Convert cloud mixing ratio into liquid and ice water path
     ! in each layer
-    do jlev = 1, nlev
-      do jcol = istartcol, iendcol
+    DO jlev = 1, nlev
+      DO jcol = istartcol, iendcol
         ! Factor to convert from gridbox-mean mass mixing ratio to
         ! in-cloud water path involves the pressure difference in
@@ -288,5 +288,5 @@
 
     if (config%iverbose >= 4) then
-      do jcol = istartcol,iendcol
+      DO jcol = istartcol,iendcol
         write(*,'(a,i0,a,f7.3,a,f7.3)') 'Profile ', jcol, ': shortwave optical depth = ', &
              &  sum(od_sw_cloud(1,:,jcol)*cloud%fraction(jcol,:)), &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_overlap.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_overlap.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_overlap.F90	(revision 5158)
@@ -94,5 +94,5 @@
     ! multiply by "factor"
     denominator = 1.0_jprb
-    do jreg = 1,nreg
+    DO jreg = 1,nreg
       op_x_frac_min(jreg) = op(jreg) &
            &  * min(frac_upper(jreg), frac_lower(jreg))
@@ -103,6 +103,6 @@
       factor = 1.0_jprb / denominator
       ! Create the random part of the overlap matrix
-      do jupper = 1,nreg
-        do jlower = 1,nreg
+      DO jupper = 1,nreg
+        DO jlower = 1,nreg
           overlap_matrix(jupper,jlower) = factor &
                &  * (frac_lower(jlower)-op_x_frac_min(jlower)) &
@@ -115,5 +115,5 @@
     
     ! Add on the maximum part of the overlap matrix
-    do jreg = 1,nreg
+    DO jreg = 1,nreg
       overlap_matrix(jreg,jreg) = overlap_matrix(jreg,jreg) &
            &  + op_x_frac_min(jreg)
@@ -368,5 +368,5 @@
 
     ! Loop through each atmospheric column
-    do jcol = istartcol, iendcol
+    DO jcol = istartcol, iendcol
       ! For this column, outer space is treated as one clear-sky
       ! region, so the fractions are assigned as such
@@ -381,5 +381,5 @@
       ! half-level starting at 1 for the top-of-atmosphere, as well
       ! as indexing each level starting at 1 for the top-most level.
-      do jlev = 1,nlev+1
+      DO jlev = 1,nlev+1
         ! Fraction of each region just below the interface
         if (jlev > nlev) then
@@ -426,6 +426,6 @@
 
         ! Convert to directional overlap matrices
-        do jupper = 1,nreg
-          do jlower = 1,nreg
+        DO jupper = 1,nreg
+          DO jlower = 1,nreg
             if (frac_lower(jlower) >= frac_threshold) then
               u_matrix(jupper,jlower,jlev,jcol) = overlap_matrix(jupper,jlower) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_pdf_sampler.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_pdf_sampler.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_pdf_sampler.F90	(revision 5158)
@@ -189,5 +189,5 @@
     real(jprb) :: wfsd, wcdf
 
-    do jsamp = 1,nsamp
+    DO jsamp = 1,nsamp
       if (mask(jsamp)) then
         ! Bilinear interpolation with bounds
@@ -237,6 +237,6 @@
     real(jprb) :: wfsd, wcdf
 
-    do jz = 1,nz
-      do jg = 1,ng
+    DO jz = 1,nz
+      DO jg = 1,ng
         if (cdf(jg, jz) > 0.0_jprb) then
           ! Bilinear interpolation with bounds
@@ -291,9 +291,9 @@
     real(jprb) :: wfsd, wcdf
 
-    do jz = 1,nz
+    DO jz = 1,nz
 
       if (mask(jz)) then
         
-        do jg = 1,ng
+        DO jg = 1,ng
           if (cdf(jg, jz) > 0.0_jprb) then
             ! Bilinear interpolation with bounds
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_random_numbers.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_random_numbers.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_random_numbers.F90	(revision 5158)
@@ -165,5 +165,5 @@
       ! populate the state
       rseed = real(abs(this%iseed),jprd)
-      do jstr = 1,this%nmaxstreams
+      DO jstr = 1,this%nmaxstreams
         ! Note that nint returns an integer of type jpib (8-byte)
         ! which may be converted to double if that is the type of
@@ -174,5 +174,5 @@
 
       ! One warmup of the C++ minstd_rand algorithm
-      do jstr = 1,this%nmaxstreams
+      DO jstr = 1,this%nmaxstreams
         this%istate(jstr) = mod(IMinstdA * this%istate(jstr), IMinstdM)
       end do
@@ -182,5 +182,5 @@
       call random_seed(size=nseed)
       allocate(iseednative(nseed))
-      do jseed = 1,nseed
+      DO jseed = 1,nseed
         iseednative(jseed) = this%iseed + jseed - 1
       end do
@@ -208,5 +208,5 @@
 
       ! C++ minstd_rand algorithm
-      do i = 1, imax
+      DO i = 1, imax
         ! The following calculation is computed entirely with 8-byte
         ! numbers (whether real or integer)
@@ -243,7 +243,7 @@
 
       ! C++ minstd_ran algorithm
-      do jblock = 1,size(randnum,2)
+      DO jblock = 1,size(randnum,2)
         ! These lines should be vectorizable
-        do i = 1, imax
+        DO i = 1, imax
           this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM)
           randnum(i,jblock) = IMinstdScale * this%istate(i)
@@ -278,8 +278,8 @@
 
       ! C++ minstd_ran algorithm
-      do jblock = 1,size(randnum,2)
+      DO jblock = 1,size(randnum,2)
         if (mask(jblock)) then
           ! These lines should be vectorizable
-          do i = 1, imax
+          DO i = 1, imax
             this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM)
             randnum(i,jblock) = IMinstdScale * this%istate(i)
@@ -290,5 +290,5 @@
     else
 
-      do jblock = 1,size(randnum,2)
+      DO jblock = 1,size(randnum,2)
         if (mask(jblock)) then
           call random_number(randnum(:,jblock))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_regions.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_regions.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_regions.F90	(revision 5158)
@@ -122,6 +122,6 @@
         ! sigma, then the 16th percentile of the lognormal is very
         ! close to exp(mu-sigma).
-        do jcol = istartcol,iendcol
-          do jlev = 1,nlev
+        DO jcol = istartcol,iendcol
+          DO jlev = 1,nlev
             if (cloud_fraction(jcol,jlev) < frac_threshold) then
               reg_fracs(1,jlev,jcol)       = 1.0_jprb
@@ -145,6 +145,6 @@
         ! cover of the denser region - see region_fractions routine
         ! below
-        do jcol = istartcol,iendcol
-          do jlev = 1,nlev
+        DO jcol = istartcol,iendcol
+          DO jlev = 1,nlev
             if (cloud_fraction(jcol,jlev) < frac_threshold) then
               reg_fracs(1,jlev,jcol)    = 1.0_jprb
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_save.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_save.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_save.F90	(revision 5158)
@@ -1168,5 +1168,5 @@
          &   dim2_name="column", dim1_name="level", &
          &   units_str="1", long_name="Ozone mass mixing ratio")
-    do jgas = 1,NMaxGases
+    DO jgas = 1,NMaxGases
       if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
         write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
@@ -1259,5 +1259,5 @@
     call gas%get(IO3, IMassMixingRatio, mixing_ratio)
     call out_file%put("o3_mmr", mixing_ratio)
-    do jgas = 1,NMaxGases
+    DO jgas = 1,NMaxGases
       if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
         write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_single_level.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_single_level.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_single_level.F90	(revision 5158)
@@ -206,5 +206,5 @@
     end if
 
-    do jcol = istartcol,iendcol
+    DO jcol = istartcol,iendcol
       this%iseed(jcol) = jcol
     end do
@@ -275,8 +275,8 @@
 
         sw_albedo_band = 0.0_jprb
-        do jband = 1,config%n_bands_sw
-          do jalbedoband = 1,nalbedoband
+        DO jband = 1,config%n_bands_sw
+          DO jalbedoband = 1,nalbedoband
             if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
-              do jcol = istartcol,iendcol
+              DO jcol = istartcol,iendcol
                 sw_albedo_band(jcol,jband) &
                     &  = sw_albedo_band(jcol,jband) & 
@@ -292,6 +292,6 @@
         if (allocated(this%sw_albedo_direct)) then
           sw_albedo_band = 0.0_jprb
-          do jband = 1,config%n_bands_sw
-            do jalbedoband = 1,nalbedoband
+          DO jband = 1,config%n_bands_sw
+            DO jalbedoband = 1,nalbedoband
               if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
                 sw_albedo_band(istartcol:iendcol,jband) &
@@ -342,8 +342,8 @@
         end if
         lw_albedo_band = 0.0_jprb
-        do jband = 1,config%n_bands_lw
-          do jalbedoband = 1,nalbedoband
+        DO jband = 1,config%n_bands_lw
+          DO jalbedoband = 1,nalbedoband
             if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
-              do jcol = istartcol,iendcol
+              DO jcol = istartcol,iendcol
                 lw_albedo_band(jcol,jband) &
                     &  = lw_albedo_band(jcol,jband) & 
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90	(revision 5158)
@@ -313,5 +313,5 @@
 
     ! Main loop over columns
-    do jcol = istartcol, iendcol
+    DO jcol = istartcol, iendcol
       ! --------------------------------------------------------
       ! Section 2: Prepare column-specific variables and arrays
@@ -325,5 +325,5 @@
       ! cloud%crop_cloud_fraction has already been called
       is_clear_sky_layer = .true.
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
           is_clear_sky_layer(jlev) = .false.
@@ -336,5 +336,5 @@
       ! In this section the reflectance, transmittance and sources
       ! are computed for each layer
-      do jlev = 1,nlev ! Start at top-of-atmosphere
+      DO jlev = 1,nlev ! Start at top-of-atmosphere
         ! --------------------------------------------------------
         ! Section 3.1: Layer-specific initialization
@@ -395,5 +395,5 @@
             ! approximate order of gas optical depth.
             ng3D = ng
-            do jg = 1, ng
+            DO jg = 1, ng
               if (od_region(jg,1) > config%max_gas_od_3D) then
                 ng3D = jg-1
@@ -488,5 +488,5 @@
               end if
 
-              do jreg = 1, nreg-1
+              DO jreg = 1, nreg-1
                 ! Compute lateral transfer rate from region jreg to
                 ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
@@ -536,5 +536,5 @@
           ! g-point, mapping from the cloud properties
           ! defined in each band.
-          do jg = 1,ng
+          DO jg = 1,ng
             ! Mapping from g-point to band
             iband = config%i_band_from_reordered_g_lw(jg)
@@ -544,5 +544,5 @@
 
             ! Loop over cloudy regions
-            do jreg = 2,nreg
+            DO jreg = 2,nreg
               ! Add scaled cloud optical depth to clear-sky value
               od_region(jg,jreg) = od_region(jg,1) &
@@ -603,5 +603,5 @@
           ! exponential, then computing the various
           ! transmission/reflectance matrices from that.
-          do jreg = 1,nregActive
+          DO jreg = 1,nregActive
             ! Write the diagonal elements of -Gamma1*z1
             Gamma_z1(1:ng3D,jreg,jreg) &
@@ -628,5 +628,5 @@
             ! non-zeros in the cloudy-sky regions even if the cloud
             ! fraction is zero
-            do jreg = nregActive+1,nreg
+            DO jreg = nregActive+1,nreg
               Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,1,1)
               Gamma_z1(1:ng3D,nreg+jreg,jreg) = Gamma_z1(1:ng3D,nreg+1,1)
@@ -656,5 +656,5 @@
           end if
 
-          do jreg = 1,nregActive-1
+          DO jreg = 1,nregActive-1
             ! Add the terms assocated with 3D transport to Gamma1*z1.
             ! First the rate from region jreg to jreg+1
@@ -775,5 +775,5 @@
           ! formulas for reflectance and transmittance, and equivalent
           ! solutions to the coupled ODEs for the sources.
-          do jreg = 2, nregActive
+          DO jreg = 2, nregActive
             call calc_reflectance_transmittance_lw(ng-ng3D, &
                  &  od_region(ng3D+1:ng,jreg), &
@@ -806,6 +806,6 @@
       ! Calculate the upwelling radiation emitted by the surface, and
       ! copy the surface albedo into total_albedo 
-      do jreg = 1,nreg
-        do jg = 1,ng
+      DO jreg = 1,nreg
+        DO jg = 1,ng
           ! region_fracs(jreg,nlev,jcol) is the fraction of each
           ! region in the lowest model level
@@ -816,5 +816,5 @@
       ! Equivalent surface values for computing clear-sky fluxes 
       if (config%do_clear) then
-        do jg = 1,ng
+        DO jg = 1,ng
           total_source_clear(jg,nlev+1) = emission(jg,jcol)
         end do
@@ -826,5 +826,5 @@
       ! Loop back up through the atmosphere computing the total albedo
       ! and the total upwelling due to emission below each level
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
         if (config%do_clear) then
           ! Use adding method for clear-sky arrays; note that there
@@ -876,5 +876,5 @@
           total_albedo_below = 0.0_jprb
           total_source_below = 0.0_jprb
-          do jreg = 1,nreg
+          DO jreg = 1,nreg
             inv_denom_scalar(:) = 1.0_jprb / (1.0_jprb &
                  &  - total_albedo(:,jreg,jreg,jlev+1)*reflectance(:,jreg,jreg,jlev))
@@ -918,6 +918,6 @@
             ! essentially diag(total_albedo) = matmul(transpose(v_matrix),
             ! diag(total_albedo_below)).
-            do jreg = 1,nreg
-              do jreg2 = 1,nreg
+            DO jreg = 1,nreg
+              DO jreg2 = 1,nreg
                 total_albedo(:,jreg,jreg,jlev) &
                      &  = total_albedo(:,jreg,jreg,jlev) &
@@ -966,5 +966,5 @@
 
       ! Final loop back down through the atmosphere to compute fluxes
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (config%do_clear) then
           ! Scalar operations for clear-sky fluxes
@@ -999,5 +999,5 @@
                &  flux_dn_above) + total_source(:,:,jlev+1)
         else
-          do jreg = 1,nreg
+          DO jreg = 1,nreg
             ! Scalar operations for all regions, requiring that
             ! reflectance, transmittance and total_albedo are diagonal
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90	(revision 5158)
@@ -328,5 +328,5 @@
 
     ! Main loop over columns
-    do jcol = istartcol, iendcol
+    DO jcol = istartcol, iendcol
       ! --------------------------------------------------------
       ! Section 2: Prepare column-specific variables and arrays
@@ -404,5 +404,5 @@
       is_clear_sky_layer = .true.
       i_cloud_top = nlev+1
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
         if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
           is_clear_sky_layer(jlev) = .false.
@@ -416,5 +416,5 @@
       ! In this section the reflectance, transmittance and sources
       ! are computed for each layer
-      do jlev = 1,nlev ! Start at top-of-atmosphere
+      DO jlev = 1,nlev ! Start at top-of-atmosphere
         ! --------------------------------------------------------
         ! Section 3.1: Layer-specific initialization
@@ -469,5 +469,5 @@
             ! approximate order of gas optical depth.
             ng3D = ng
-            do jg = 1, ng
+            DO jg = 1, ng
               if (od_region(jg,1) > config%max_gas_od_3D) then
                 ng3D = jg-1
@@ -553,5 +553,5 @@
               end if
 
-              do jreg = 1, nreg-1
+              DO jreg = 1, nreg-1
                 ! Compute lateral transfer rates from region jreg to
                 ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
@@ -615,5 +615,5 @@
           ! g-point, mapping from the cloud properties
           ! defined in each band.
-          do jg = 1,ng
+          DO jg = 1,ng
             ! Mapping from g-point to band
             iband = config%i_band_from_reordered_g_sw(jg)
@@ -629,5 +629,5 @@
 
             ! Loop over cloudy regions
-            do jreg = 2,nreg
+            DO jreg = 2,nreg
               scat_od_cloud = od_cloud(iband,jlev,jcol) &
                    &  * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
@@ -679,5 +679,5 @@
           ! exponential, then computing the various
           ! transmission/reflectance matrices from that.
-          do jreg = 1,nregactive
+          DO jreg = 1,nregactive
             ! Write the diagonal elements of -Gamma1*z1
             Gamma_z1(1:ng3D,jreg,jreg) &
@@ -701,5 +701,5 @@
           end do
 
-          do jreg = 1,nregactive-1
+          DO jreg = 1,nregactive-1
             ! Write the elements of -Gamma1*z1 concerned with 3D
             ! transport
@@ -820,5 +820,5 @@
           ! Compute reflectance, transmittance and associated terms
           ! for each cloudy region, using the Meador-Weaver formulas
-          do jreg = 2, nregactive
+          DO jreg = 2, nregactive
             call calc_reflectance_transmittance_sw(ng-ng3D, &
                  &  mu0, &
@@ -853,6 +853,6 @@
       ! beam incident on the surface, and copy the surface albedo
       ! into total_albedo
-      do jreg = 1,nreg
-        do jg = 1,ng
+      DO jreg = 1,nreg
+        DO jg = 1,ng
           total_albedo(jg,jreg,jreg,nlev+1) = albedo_diffuse(jg,jcol)
           total_albedo_direct(jg,jreg,jreg,nlev+1) &
@@ -875,5 +875,5 @@
       ! Work back up through the atmosphere computing the total albedo
       ! of the atmosphere below that point using the adding method
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
 
         ! --------------------------------------------------------
@@ -995,6 +995,6 @@
           ! radiation:
           total_albedo(:,:,:,jlev) = 0.0_jprb
-          do jreg = 1,nreg    ! Target layer (jlev-1)
-            do jreg2 = 1,nreg ! Current layer (jlev)
+          DO jreg = 1,nreg    ! Target layer (jlev-1)
+            DO jreg2 = 1,nreg ! Current layer (jlev)
               total_albedo(:,jreg,jreg,jlev) = total_albedo(:,jreg,jreg,jlev) &
                    &  + sum(total_albedo_below(:,:,jreg2),2) &
@@ -1004,6 +1004,6 @@
           ! ...then direct radiation:
           total_albedo_direct(:,:,:,jlev) = 0.0_jprb
-          do jreg = 1,nreg    ! Target layer (jlev-1)
-            do jreg2 = 1,nreg ! Current layer (jlev)
+          DO jreg = 1,nreg    ! Target layer (jlev-1)
+            DO jreg2 = 1,nreg ! Current layer (jlev)
               total_albedo_direct(:,jreg,jreg,jlev) = total_albedo_direct(:,jreg,jreg,jlev) &
                    &  + sum(total_albedo_below_direct(:,:,jreg2),2) &
@@ -1031,5 +1031,5 @@
           ! First diffuse radiation:
           albedo_part = total_albedo_below
-          do jreg = 1,nreg
+          DO jreg = 1,nreg
             albedo_part(:,jreg,jreg) = 0.0_jprb
           end do
@@ -1040,5 +1040,5 @@
           ! ...then direct radiation:
           albedo_part = total_albedo_below_direct
-          do jreg = 1,nreg
+          DO jreg = 1,nreg
             albedo_part(:,jreg,jreg) = 0.0_jprb
           end do
@@ -1059,6 +1059,6 @@
             ! essentially diag(total_albedo) += matmul(transpose(v_matrix),
             ! diag(total_albedo_below)).
-            do jreg = 1,nreg
-              do jreg2 = 1,nreg
+            DO jreg = 1,nreg
+              DO jreg2 = 1,nreg
                 total_albedo(:,jreg,jreg,jlev) &
                      &  = total_albedo(:,jreg,jreg,jlev) &
@@ -1075,5 +1075,5 @@
             ! "Explicit entrapment"
 
-            do jreg2 = 1,nreg
+            DO jreg2 = 1,nreg
               ! Loop through each region in the lower layer. For one
               ! of the regions in the lower layer, we are imagining it
@@ -1110,5 +1110,5 @@
                      &  / max(config%cloud_fraction_threshold, region_fracs(jreg2,jlev,jcol))
 
-                do jreg = 1, nreg-1
+                DO jreg = 1, nreg-1
                   ! Compute lateral transfer rates from region jreg to
                   ! jreg+1 as before, but without length scale which
@@ -1139,5 +1139,5 @@
               inv_effective_size = min(cloud%inv_cloud_effective_size(jcol,jlev-1), &
                    &                   1.0_jprb/config%min_cloud_effective_size)
-              do jreg = 1,nreg-1
+              DO jreg = 1,nreg-1
                 ! Diffuse transport down and up with one random
                 ! scattering event
@@ -1167,5 +1167,5 @@
               ! trival limit, so we cap the maximum input to expm by
               ! scaling down if necessary
-              do jg = 1,ng
+              DO jg = 1,ng
                 max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2))
                 if (max_entr > config%max_cloud_od) then
@@ -1197,6 +1197,6 @@
 #ifndef EXPLICIT_EDGE_ENTRAPMENT
               ! Scale to get the contribution to the diffuse albedo
-              do jreg3 = 1,nreg
-                do jreg = 1,nreg
+              DO jreg3 = 1,nreg
+                DO jreg = 1,nreg
                   albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) &
                        &  * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg2)
@@ -1210,11 +1210,11 @@
               albedo_part = 0.0_jprb
               ! Scale to get the contribution to the diffuse albedo
-              do jreg3 = 1,nreg     ! TO upper region
-                do jreg = 1,nreg    ! FROM upper region
+              DO jreg3 = 1,nreg     ! TO upper region
+                DO jreg = 1,nreg    ! FROM upper region
                   transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & 
                        &  * cloud%overlap_param(jcol,jlev-1) &
                        &  * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev,jcol)) &
                        &  / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol))
-                  do jreg4 = 1,nreg ! VIA first lower region (jreg2 is second lower region)
+                  DO jreg4 = 1,nreg ! VIA first lower region (jreg2 is second lower region)
                     if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then
                       albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) &
@@ -1236,5 +1236,5 @@
               ! Now do the same for the direct albedo
               entrapment = 0.0_jprb
-              do jreg = 1,nreg-1
+              DO jreg = 1,nreg-1
                 ! Direct transport down and diffuse up with one random
                 ! scattering event
@@ -1264,5 +1264,5 @@
               ! trival limit, so we cap the maximum input to expm by
               ! scaling down if necessary
-              do jg = 1,ng
+              DO jg = 1,ng
                 max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2))
                 if (max_entr > config%max_cloud_od) then
@@ -1289,6 +1289,6 @@
 
 #ifndef EXPLICIT_EDGE_ENTRAPMENT
-              do jreg3 = 1,nreg
-                do jreg = 1,nreg
+              DO jreg3 = 1,nreg
+                DO jreg = 1,nreg
                   albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) &
                        &  * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg2)
@@ -1298,11 +1298,11 @@
               entrapment = albedo_part
               albedo_part = 0.0_jprb
-              do jreg3 = 1,nreg
-                do jreg = 1,nreg
+              DO jreg3 = 1,nreg
+                DO jreg = 1,nreg
                   transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & 
                        &  * cloud%overlap_param(jcol,jlev-1) &
                        &  * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev-1,jcol)) &
                        &  / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol))
-                  do jreg4 = 1,nreg
+                  DO jreg4 = 1,nreg
                     if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then
                      albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) &
@@ -1343,6 +1343,6 @@
           end if
 
-          do jreg = 1,nreg          ! Target layer (jlev-1)
-            do jreg2 = 1,nregactive ! Current layer (jlev)
+          DO jreg = 1,nreg          ! Target layer (jlev-1)
+            DO jreg2 = 1,nregactive ! Current layer (jlev)
               x_direct_above(:,jreg) = x_direct_above(:,jreg) &
                    &  + x_direct(:,jreg2) * v_matrix(jreg2,jreg,jlev,jcol)
@@ -1368,5 +1368,5 @@
       ! Direct downwelling flux (into a plane perpendicular to the
       ! sun) entering the top of each region in the top-most layer
-      do jreg = 1,nreg
+      DO jreg = 1,nreg
         direct_dn_below(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
       end do
@@ -1421,5 +1421,5 @@
 
       ! Final loop back down through the atmosphere to compute fluxes
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
 
 #ifdef PRINT_ENTRAPMENT_DATA
@@ -1675,5 +1675,5 @@
          &                             + tan_diffuse_angle_3d*tan_diffuse_angle_3d) * 0.5_jprb
 
-    do jreg = istartreg,iendreg
+    DO jreg = istartreg,iendreg
       ! Geometric series enhancement due to multiple scattering: the
       ! amplitude enhancement is equal to the limit of
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90	(revision 5158)
@@ -201,5 +201,5 @@
     else
       find_wavenumber = 1
-      do while (wavenumber > this%wavenumber2(find_wavenumber) &
+      DO while (wavenumber > this%wavenumber2(find_wavenumber) &
            &    .and. find_wavenumber < this%nwav)
         find_wavenumber = find_wavenumber + 1
@@ -284,7 +284,7 @@
       end if
 
-      do jband = 1,this%nband
+      DO jband = 1,this%nband
         weight = 0.0_jprb
-        do jwav = 1,nwav
+        DO jwav = 1,nwav
           ! Work out wavenumber range for which this cloud wavenumber
           ! will be applicable
@@ -321,5 +321,5 @@
             ! Find interpolating points
             iwav = 2
-            do while (wavenumber(iwav) < this%wavenumber2_band(jband))
+            DO while (wavenumber(iwav) < this%wavenumber2_band(jband))
               iwav = iwav+1
             end do
@@ -358,5 +358,5 @@
       mapping = 0.0_jprb
       ! Loop over wavenumbers representing cloud
-      do jwav = 1,nwav
+      DO jwav = 1,nwav
         ! Clear the weights. The weight says for one wavenumber in the
         ! cloud file, what is its fractional contribution to each of
@@ -400,5 +400,5 @@
                  &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
             if (isd1-isd0 > 1) then
-              do isd = isd0+1,isd1-1
+              DO isd = isd0+1,isd1-1
                 ! Intermediate trapezia
                 weight(isd) = 0.5_jprb * (this%wavenumber1(isd)+this%wavenumber2(isd) &
@@ -447,5 +447,5 @@
                  &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
             if (isd2-isd1 > 1) then
-              do isd = isd1+1,isd2-1
+              DO isd = isd1+1,isd2-1
                 ! Intermediate trapezia
                 weight(isd) = weight(isd) + 0.5_jprb * (2.0_jprb*wavenum2 &
@@ -468,5 +468,5 @@
         weight = weight * planck_weight
 
-        do jg = 1,this%ng
+        DO jg = 1,this%ng
           mapping(jg, jwav) = sum(weight * this%gpoint_fraction(:,jg))
         end do
@@ -478,5 +478,5 @@
 
       ! Normalize mapping matrix
-      do jg = 1,this%ng
+      DO jg = 1,this%ng
         mapping(jg,:) = mapping(jg,:) * (1.0_jprb/sum(mapping(jg,:)))
       end do
@@ -589,5 +589,5 @@
     ! Check wavelength is monotonically increasing
     if (ninterval > 2) then
-      do jint = 2,ninterval-1
+      DO jint = 2,ninterval-1
         if (wavelength_bound(jint) <= wavelength_bound(jint-1)) then
           write(nulerr, '(a)') '*** Error: wavelength bounds must be monotonically increasing'
@@ -609,6 +609,6 @@
       end if
 
-      do jband = 1,this%nband
-        do jint = 1,ninterval
+      DO jband = 1,this%nband
+        DO jint = 1,ninterval
           if (jint == 1) then
             ! First input interval in wavelength space: lower
@@ -692,5 +692,5 @@
       
       ! All bounded intervals
-      do jint = 2,ninterval-1
+      DO jint = 2,ninterval-1
         wavenumber1_bound = 0.01_jprb / wavelength_bound(jint)
         wavenumber2_bound = 0.01_jprb / wavelength_bound(jint-1)
@@ -710,6 +710,6 @@
       end if
 
-      do jg = 1,this%ng
-        do jin = 1,ninput
+      DO jg = 1,this%ng
+        DO jin = 1,ninput
           mapping(jin,jg) = sum(this%gpoint_fraction(:,jg) * planck, &
                &                 mask=(i_input==jin))
@@ -723,7 +723,7 @@
 
       ! Loop through all intervals
-      do jint = 1,ninterval
+      DO jint = 1,ninterval
         ! Loop through the wavenumbers for gpoint_fraction
-        do jwav = 1,this%nwav
+        DO jwav = 1,this%nwav
           if (jint == 1) then
             ! First input interval in wavelength space: lower
@@ -758,5 +758,5 @@
       end do
       if (use_fluxes_local) then
-        do jg = 1,this%ng
+        DO jg = 1,this%ng
           mapping(:,jg) = mapping(:,jg) / sum(this%gpoint_fraction(:,jg) * planck)
         end do
@@ -769,5 +769,5 @@
     if (.not. use_fluxes_local) then
       ! Normalize mapping matrix
-      do jg = 1,size(mapping,dim=2)
+      DO jg = 1,size(mapping,dim=2)
         mapping(:,jg) = mapping(:,jg) * (1.0_jprb/sum(mapping(:,jg)))
       end do
@@ -828,5 +828,5 @@
     is_band_unassigned = .true.
 
-    do jint = 1,ninterval
+    DO jint = 1,ninterval
       inext = minloc(wavelength1, dim=1, mask=is_band_unassigned)
       if (jint > 1) then
@@ -873,8 +873,8 @@
       write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' bands (wavenumber ranges in cm-1)'
       if (nout <= 40) then
-        do jout = 1,nout
+        DO jout = 1,nout
           write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
                &                        nint(this%wavenumber2_band(jout)), ':'
-          do jin = 1,nin
+          DO jin = 1,nin
             write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
           end do
@@ -882,8 +882,8 @@
         end do
       else
-        do jout = 1,30
+        DO jout = 1,30
           write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
                &                        nint(this%wavenumber2_band(jout)), ':'
-          do jin = 1,nin
+          DO jin = 1,nin
             write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
           end do
@@ -893,5 +893,5 @@
         write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(nout)), ' to', &
              &                        nint(this%wavenumber2_band(nout)), ':'
-        do jin = 1,nin
+        DO jin = 1,nin
           write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
         end do
@@ -901,7 +901,7 @@
       write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' g-points'
       if (nout <= 40) then
-        do jout = 1,nout
+        DO jout = 1,nout
           write(nulout,'(i3,a)',advance='no') jout, ':'
-          do jin = 1,nin
+          DO jin = 1,nin
             write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
           end do
@@ -909,7 +909,7 @@
         end do
       else
-        do jout = 1,30
+        DO jout = 1,30
           write(nulout,'(i3,a)',advance='no') jout, ':'
-          do jin = 1,nin
+          DO jin = 1,nin
             write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
           end do
@@ -918,5 +918,5 @@
         write(nulout,'(a)') '  ...'
         write(nulout,'(i3,a)',advance='no') nout, ':'
-        do jin = 1,nin
+        DO jin = 1,nin
           write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
         end do
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_thermodynamics.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_thermodynamics.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_thermodynamics.F90	(revision 5158)
@@ -143,6 +143,6 @@
     end if
 
-    do jlev = 1,nlev
-       do jcol = istartcol,iendcol
+    DO jlev = 1,nlev
+       DO jcol = istartcol,iendcol
           pressure = 0.5 * (this%pressure_hl(jcol,jlev)+this%pressure_hl(jcol,jlev+1))
           temperature = 0.5 * (this%temperature_hl(jcol,jlev)+this%temperature_hl(jcol,jlev+1))
@@ -265,5 +265,5 @@
       ! be half the separation between the half-levels at the edge of
       ! the two adjacent layers
-      do jlev = 2,nlev-1
+      DO jlev = 2,nlev-1
         layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
              &                    * log(pressure_hl(i1:i2,jlev+2)/pressure_hl(i1:i2,jlev))
@@ -276,5 +276,5 @@
       ! don't take the logarithm of the last pressure in each column.
 
-      do jlev = 1,nlev-2
+      DO jlev = 1,nlev-2
         layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
              &                    * log(pressure_hl(i1:i2,jlev)/pressure_hl(i1:i2,jlev+2))
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90	(revision 5158)
@@ -202,5 +202,5 @@
 
     ! Main loop over columns
-    do jcol = istartcol, iendcol
+    DO jcol = istartcol, iendcol
       ! --------------------------------------------------------
       ! Section 2: Prepare column-specific variables and arrays
@@ -211,5 +211,5 @@
       is_clear_sky_layer = .true.
       i_cloud_top = nlev+1
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
           is_clear_sky_layer(jlev) = .false.
@@ -268,5 +268,5 @@
         ! Save the spectral fluxes if required
         if (config%do_save_spectral_flux) then
-          do jlev = 1,nlev+1
+          DO jlev = 1,nlev+1
             call indexed_sum(flux_up_clear(:,jlev), &
                  &           config%i_spec_from_reordered_g_lw, &
@@ -294,5 +294,5 @@
       end if
 
-      do jlev = i_cloud_top,nlev ! Start at cloud top and work down
+      DO jlev = i_cloud_top,nlev ! Start at cloud top and work down
 
         ! Copy over clear-sky properties
@@ -308,5 +308,5 @@
           source_dn(:,2:,jlev)     = 0.0_jprb
         else
-          do jreg = 2,nreg
+          DO jreg = 2,nreg
             ! Cloudy sky
             ! Add scaled cloud optical depth to clear-sky value
@@ -358,5 +358,5 @@
           end do
           ! Emission is scaled by the size of each region
-          do jreg = 1,nregions
+          DO jreg = 1,nregions
             source_up(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_up(:,jreg,jlev)
             source_dn(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_dn(:,jreg,jlev)
@@ -375,6 +375,6 @@
       ! Calculate the upwelling radiation emitted by the surface, and
       ! copy the surface albedo into total_albedo 
-      do jreg = 1,nregions
-        do jg = 1,ng
+      DO jreg = 1,nregions
+        DO jg = 1,ng
           ! region_fracs(jreg,nlev,jcol) is the fraction of each region in the
           ! lowest model level
@@ -387,5 +387,5 @@
       ! atmosphere and the total upwelling due to emission below each
       ! level below using the adding method
-      do jlev = nlev,i_cloud_top,-1
+      DO jlev = nlev,i_cloud_top,-1
 
         total_albedo_below        = 0.0_jprb
@@ -394,5 +394,5 @@
           total_albedo_below = 0.0_jprb
           total_source_below = 0.0_jprb
-          do jg = 1,ng
+          DO jg = 1,ng
             inv_denom(jg,1) = 1.0_jprb &
                  &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev))
@@ -428,6 +428,6 @@
           ! the operation we perform is essentially diag(total_albedo)
           ! = matmul(transpose(v_matrix), diag(total_albedo_below)).
-          do jreg = 1,nregions
-            do jreg2 = 1,nregions
+          DO jreg = 1,nregions
+            DO jreg2 = 1,nregions
               total_albedo(:,jreg,jlev) &
                    &  = total_albedo(:,jreg,jlev) &
@@ -445,5 +445,5 @@
       ! Section 6: Copy over downwelling fluxes above cloud top
       ! --------------------------------------------------------
-      do jlev = 1,i_cloud_top
+      DO jlev = 1,i_cloud_top
         if (config%do_clear) then
           ! Clear-sky fluxes have already been averaged: use these
@@ -476,5 +476,5 @@
              &           flux%lw_up_band(:,jcol,i_cloud_top))
       end if
-      do jlev = i_cloud_top-1,1,-1
+      DO jlev = i_cloud_top-1,1,-1
         flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
         flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
@@ -494,13 +494,13 @@
       ! scattering layer, using overlap matrix to translate to the
       ! regions of the first layer of cloud
-      do jreg = 1,nregions
+      DO jreg = 1,nregions
         flux_dn(:,jreg)  = v_matrix(jreg,1,i_cloud_top,jcol)*flux_dn_clear(:,i_cloud_top)
       end do
 
       ! Final loop back down through the atmosphere to compute fluxes
-      do jlev = i_cloud_top,nlev
+      DO jlev = i_cloud_top,nlev
 
         if (is_clear_sky_layer(jlev)) then
-          do jg = 1,ng
+          DO jg = 1,ng
             flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) &
                  &  + reflectance(jg,1,jlev)*total_source(jg,1,jlev+1) + source_dn(jg,1,jlev) ) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90	(revision 5158)
@@ -127,5 +127,5 @@
     ! entire clear-sky atmosphere at once
     real(jprb), dimension(config%n_g_sw, nlev) &
-    	 &  :: reflectance_clear, transmittance_clear, &
+         &  :: reflectance_clear, transmittance_clear, &
          &     ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear
 
@@ -198,5 +198,5 @@
 
     ! Main loop over columns
-    do jcol = istartcol, iendcol
+    DO jcol = istartcol, iendcol
       ! --------------------------------------------------------
       ! Section 2: Prepare column-specific variables and arrays
@@ -251,5 +251,5 @@
       ! cloud%crop_cloud_fraction has already been called
       is_clear_sky_layer = .true.
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
           is_clear_sky_layer(jlev) = .false.
@@ -271,8 +271,8 @@
 
       ! Cloudy layers
-      do jlev = 1,nlev ! Start at top-of-atmosphere
+      DO jlev = 1,nlev ! Start at top-of-atmosphere
         if (.not. is_clear_sky_layer(jlev)) then
-          do jreg = 2,nregions
-            do jg = 1,ng
+          DO jreg = 2,nregions
+            DO jg = 1,ng
               ! Mapping from g-point to band
               iband = config%i_band_from_reordered_g_sw(jg)
@@ -315,5 +315,5 @@
 
       ! Copy surface albedos in clear-sky region
-      do jg = 1,ng
+      DO jg = 1,ng
         total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
         total_albedo_direct(jg,1,nlev+1) &
@@ -324,5 +324,5 @@
       ! underneath
       if (.not. is_clear_sky_layer(nlev)) then
-        do jreg = 2,nregions
+        DO jreg = 2,nregions
           total_albedo(:,jreg,nlev+1)        = total_albedo(:,1,nlev+1)
           total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1)
@@ -337,5 +337,5 @@
       ! Work up from the surface computing the total albedo of the
       ! atmosphere below that point using the adding method
-      do jlev = nlev,1,-1
+      DO jlev = nlev,1,-1
 
         total_albedo_below        = 0.0_jprb
@@ -346,5 +346,5 @@
           ! "below" quantities since with no cloud overlap to worry
           ! about, these are the same
-          do jg = 1,ng
+          DO jg = 1,ng
             inv_denom(jg,1) = 1.0_jprb &
                  &  / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev))
@@ -361,5 +361,5 @@
 
         ! All-sky fluxes: first the clear region
-        do jg = 1,ng
+        DO jg = 1,ng
           inv_denom(jg,1) = 1.0_jprb &
                &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev))
@@ -375,6 +375,6 @@
         ! Then the cloudy regions if any cloud is present in this layer
         if (.not. is_clear_sky_layer(jlev)) then
-          do jreg = 2,nregions
-            do jg = 1,ng
+          DO jreg = 2,nregions
+            DO jg = 1,ng
               inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb &
                    &  - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev))
@@ -400,6 +400,6 @@
           ! the operation we perform is essentially diag(total_albedo)
           ! = matmul(transpose(v_matrix)), diag(total_albedo_below)).
-          do jreg = 1,nregions
-            do jreg2 = 1,nregions
+          DO jreg = 1,nregions
+            DO jreg2 = 1,nregions
               total_albedo(:,jreg,jlev) &
                    &  = total_albedo(:,jreg,jlev) &
@@ -426,5 +426,5 @@
       ! Direct downwelling flux (into a plane perpendicular to the
       ! sun) entering the top of each region in the top-most layer
-      do jreg = 1,nregions
+      DO jreg = 1,nregions
         direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
       end do
@@ -495,7 +495,7 @@
 
       ! Final loop back down through the atmosphere to compute fluxes
-      do jlev = 1,nlev
+      DO jlev = 1,nlev
         if (config%do_clear) then
-          do jg = 1,ng
+          DO jg = 1,ng
             flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
                &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1)*reflectance_clear(jg,jlev) &
@@ -509,5 +509,5 @@
 
         ! All-sky fluxes: first the clear region...
-        do jg = 1,ng
+        DO jg = 1,ng
           flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
                &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1)*reflectance_clear(jg,jlev) &
@@ -525,6 +525,6 @@
           direct_dn(:,2:nregions)= 0.0_jprb
         else
-          do jreg = 2,nregions
-            do jg = 1,ng
+          DO jreg = 2,nregions
+            DO jg = 1,ng
               flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) &
                    &  * (trans_dir_dir(jg,jreg,jlev)*total_albedo_direct(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev) &
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_two_stream.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_two_stream.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_two_stream.F90	(revision 5158)
@@ -72,5 +72,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jg = 1, ng
+    DO jg = 1, ng
       ! Fu et al. (1997), Eq 2.9 and 2.10:
       !      gamma1(jg) = LwDiffusivity * (1.0_jprb - 0.5_jprb*ssa(jg) &
@@ -122,5 +122,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jg = 1, ng
+    DO jg = 1, ng
       !      gamma1(jg) = 2.0_jprb  - ssa(jg) * (1.25_jprb + 0.75_jprb*g(jg))
       !      gamma2(jg) = 0.75_jprb *(ssa(jg) * (1.0_jprb - g(jg)))
@@ -195,5 +195,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jg = 1, ng
+    DO jg = 1, ng
       if (od(jg) > 1.0e-3_jprd) then
         k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
@@ -291,5 +291,5 @@
 #endif
 
-    do jg = 1, ng
+    DO jg = 1, ng
       factor = (LwDiffusivityWP * 0.5_jprb) * ssa(jg)
       gamma1 = LwDiffusivityWP - factor*(1.0_jprb + asymmetry(jg))
@@ -379,5 +379,5 @@
 #endif
 
-    do jg = 1, ng
+    DO jg = 1, ng
       ! Compute upward and downward emission assuming the Planck
       ! function to vary linearly with optical depth within the layer
@@ -474,5 +474,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jg = 1, ng
+    DO jg = 1, ng
 
         gamma4 = 1.0_jprd - gamma3(jg)
@@ -624,5 +624,5 @@
     trans_dir_dir = exp(trans_dir_dir)
 
-    do jg = 1, ng
+    DO jg = 1, ng
 
       ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to
@@ -652,5 +652,5 @@
     exponential = exp(-k_exponent*od)
 
-    do jg = 1, ng
+    DO jg = 1, ng
       k_mu0 = k_exponent(jg)*mu0
       one_minus_kmu0_sqr = 1.0_jprb - k_mu0*k_mu0
@@ -698,5 +698,5 @@
 #else
     ! GPU-capable and vector-optimized version for ICON
-    do jg = 1, ng
+    DO jg = 1, ng
 
       trans_dir_dir(jg) = max(-max(od(jg) * (1.0_jprb/mu0),0.0_jprb),-1000.0_jprb)
@@ -810,5 +810,5 @@
 ! Added for DWD (2020)
 !NEC$ shortloop
-    do jg = 1, ng
+    DO jg = 1, ng
       ! Note that if the minimum value is reduced (e.g. to 1.0e-24)
       ! then noise starts to appear as a function of solar zenith
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/easy_netcdf.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/easy_netcdf.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/easy_netcdf.F90	(revision 5158)
@@ -392,5 +392,5 @@
 
     ndimlens(:) = 0
-    do j = 1,ndims
+    DO j = 1,ndims
       istatus = nf90_inquire_dimension(this%ncid, idimids(j), len=ndimlens(j))
       if (istatus /= NF90_NOERR) then
@@ -404,5 +404,5 @@
     if (present(ntotal)) then
       ntotal = 1
-      do j = 1, ndims
+      DO j = 1, ndims
         ntotal = ntotal * ndimlens(j)
       end do
@@ -613,5 +613,5 @@
     ! Compute number of elements of the variable in the file
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
     end do
@@ -658,5 +658,5 @@
     ! Compute number of elements of the variable in the file
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
     end do
@@ -707,5 +707,5 @@
     ! i.e. excluding the slowest varying dimension, indexed by "index"
     ntotal = 1
-    do j = 1, ndims-1
+    DO j = 1, ndims-1
       ntotal = ntotal * ndimlens(j)
     end do
@@ -761,5 +761,5 @@
     ! Ensure variable has only one dimension in the file
     n = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       n = n * ndimlens(j)
       if (j > 1 .and. ndimlens(j) > 1) then
@@ -819,5 +819,5 @@
     ! Ensure variable has only one dimension in the file
     n = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       n = n * ndimlens(j)
       if (j > 1 .and. ndimlens(j) > 1) then
@@ -878,5 +878,5 @@
     ! Ensure variable has only one dimension in the file
     n = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       n = n * ndimlens(j)
       if (j > 1 .and. ndimlens(j) > 1) then
@@ -938,5 +938,5 @@
     ! Ensure variable has only one dimension aside from the last one
     n = 1
-    do j = 1, ndims-1
+    DO j = 1, ndims-1
       n = n * ndimlens(j)
       if (j > 1 .and. ndimlens(j) > 1) then
@@ -1021,5 +1021,5 @@
     ! dimensions
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
       if (j > 2 .and. ndimlens(j) > 1) then
@@ -1133,5 +1133,5 @@
     ! dimensions
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
       if (j > 2 .and. ndimlens(j) > 1) then
@@ -1202,5 +1202,5 @@
       vcount(1:2) = [ndimlen1,1]
       
-      do j = 1,ndimlen2
+      DO j = 1,ndimlen2
         vstart(2) = j
         istatus = nf90_get_var(this%ncid, ivarid, matrix(:,j), start=vstart, count=vcount)
@@ -1252,5 +1252,5 @@
     ! dimensions aside from the last one
     ntotal = 1
-    do j = 1, ndims-1
+    DO j = 1, ndims-1
       ntotal = ntotal * ndimlens(j)
       if (j > 2 .and. ndimlens(j) > 1) then
@@ -1376,5 +1376,5 @@
     ! dimensions
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
       if (j > 3 .and. ndimlens(j) > 1) then
@@ -1512,5 +1512,5 @@
     ! dimensions aside from the last one
     ntotal = 1
-    do j = 1, ndims-1
+    DO j = 1, ndims-1
       ntotal = ntotal * ndimlens(j)
       if (j > 3 .and. ndimlens(j) > 1) then
@@ -1654,5 +1654,5 @@
     ! dimensions
     ntotal = 1
-    do j = 1, ndims
+    DO j = 1, ndims
       ntotal = ntotal * ndimlens(j)
       if (j > 4 .and. ndimlens(j) > 1) then
@@ -1804,5 +1804,5 @@
 
     ! Pad with blanks since nf90_get_att does not do this
-    do j = i_attr_len+1,len(attr_str)
+    DO j = i_attr_len+1,len(attr_str)
       attr_str(j:j) = ' '
     end do
@@ -1878,5 +1878,5 @@
 
     ! Pad with blanks since nf90_get_att does not do this
-    do j = i_attr_len+1,len(attr_str)
+    DO j = i_attr_len+1,len(attr_str)
       attr_str(j:j) = ' '
     end do
@@ -2705,5 +2705,5 @@
     end if
 
-    do jdim = 1,ndims
+    DO jdim = 1,ndims
       istatus = nf90_inquire_dimension(infile%ncid, idimids(jdim), &
            &  name=dimname, len=dimlen)
@@ -2771,5 +2771,5 @@
 
     ! Map dimension IDs
-    do jdim = 1,ndims
+    DO jdim = 1,ndims
       istatus = nf90_inquire_dimension(infile%ncid, idimids_in(jdim), name=dim_name)
       if (istatus /= NF90_NOERR) then
@@ -2801,5 +2801,5 @@
 
     ! Copy attributes
-    do jattr = 1,nattr
+    DO jattr = 1,nattr
       istatus = nf90_inq_attname(infile%ncid, ivarid_in, jattr, attr_name)
       if (istatus /= NF90_NOERR) then
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/print_matrix.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/print_matrix.F90	(revision 5154)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/print_matrix.F90	(revision 5158)
@@ -16,5 +16,5 @@
       write(unit_local,'(a,a,$)') name, '=['
     end if
-    do i = 1,size(vec,1)
+    DO i = 1,size(vec,1)
        write(unit_local,'(f16.8,$)') vec(i)
     end do
@@ -43,6 +43,6 @@
       write(unit_local,'(a,a,$)') name, '=['
     end if
-    do i = 1,size(mat,1)
-       do j = 1,size(mat,2)
+    DO i = 1,size(mat,1)
+       DO j = 1,size(mat,2)
           write(unit_local,'(f16.8,$)') mat(i,j)
        end do
@@ -61,7 +61,7 @@
     integer    :: i, j, k
     write(6,'(a,a)') name, '='
-    do k = 1,size(mat,3)
-      do i = 1,size(mat,1)
-        do j = 1,size(mat,2)
+    DO k = 1,size(mat,3)
+      DO i = 1,size(mat,1)
+        DO j = 1,size(mat,2)
              write(6,'(f16.8,$)') mat(i,j,k)
           end do
