Index: /LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2691)
+++ /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2692)
@@ -1090,5 +1090,5 @@
     ! en imposant la valeur de igout.
     !======================================================================d
-    if (prt_level.ge.1) then
+    IF (prt_level.ge.1) THEN
        igout=klon/2+1/klon
        write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
@@ -1101,26 +1101,26 @@
 
        write(lunout,*) 'paprs, play, phi, u, v, t'
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
                u(igout,k),v(igout,k),t(igout,k)
-       enddo
+       ENDDO
        write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
-       enddo
-    endif
+       ENDDO
+    ENDIF
 
     !======================================================================
 
-    if (first) then 
+    IF (first) THEN 
        !CR:nvelles variables convection/poches froides
 
        print*, '================================================='
        print*, 'Allocation des variables locales et sauvegardees'
-       call phys_local_var_init
+       CALL phys_local_var_init
        !
        pasphys=pdtphys
        !     appel a la lecture du run.def physique
-       call conf_phys(ok_journe, ok_mensuel, &
+       CALL conf_phys(ok_journe, ok_mensuel, &
             ok_instan, ok_hf, &
             ok_LES, &
@@ -1136,14 +1136,14 @@
             read_climoz, &
             alp_offset)
-       call phys_state_var_init(read_climoz)
-       call phys_output_var_init
+       CALL phys_state_var_init(read_climoz)
+       CALL phys_output_var_init
        print*, '================================================='
        !
        !CR: check sur le nb de traceurs de l eau
-       if ((iflag_ice_thermo.gt.0).and.(nqo==2)) then
+       IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
           WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
                '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
           STOP
-       endif
+       ENDIF
 
        dnwd0=0.0
@@ -1158,5 +1158,5 @@
        first=.false.
 
-    endif  ! first
+    ENDIF  ! first
 
     !ym => necessaire pour iflag_con != 2    
@@ -1180,6 +1180,6 @@
        DO i=1,klon
           zero_v(i)=0.
-       END DO
-    END IF
+       ENDDO
+    ENDIF
 
     IF (debut) THEN
@@ -1195,5 +1195,5 @@
     ENDIF
 
-    if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
+    IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
 
 
@@ -1239,5 +1239,5 @@
        ELSE 
           config_inca='none' ! default
-       END IF
+       ENDIF
 
        IF (aerosol_couple .AND. (config_inca /= "aero" &
@@ -1278,7 +1278,7 @@
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-       if (iflag_pbl>1) then
+       IF (iflag_pbl>1) THEN
           PRINT*, "Using method MELLOR&YAMADA" 
-       endif
+       ENDIF
 
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1298,5 +1298,5 @@
           abort_message='nbre de pas de temps physique n est pas multiple ' &
                // 'de nbapp_rad'
-          call abort_physic(modname,abort_message,1)
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1320,5 +1320,4 @@
 
 
-
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !
@@ -1343,5 +1342,5 @@
                klon
           abort_message='nlon et klon ne sont pas coherents'
-          call abort_physic(modname,abort_message,1)
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
        IF (nlev .NE. klev) THEN
@@ -1349,5 +1348,5 @@
                klev
           abort_message='nlev et klev ne sont pas coherents'
-          call abort_physic(modname,abort_message,1)
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
        !
@@ -1356,5 +1355,5 @@
           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
           abort_message='Nbre d appels au rayonnement insuffisant'
-          call abort_physic(modname,abort_message,1)
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
        WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
@@ -1387,5 +1386,5 @@
           !CR:04.12.07: initialisations poches froides
           ! Controle de ALE et ALP pour la fermeture convective (jyg)
-          if (iflag_wake>=1) then
+          IF (iflag_wake>=1) THEN
              CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
                   ,alp_bl_prescr, ale_bl_prescr)
@@ -1407,5 +1406,5 @@
              d_s_wk(:) = 0.
              d_dens_wk(:) = 0.
-          endif
+          ENDIF
 
           !        do i = 1,klon
@@ -1419,5 +1418,5 @@
           OPEN(98,file='npCFMIP_param.data',status='old', &
                form='formatted',iostat=iostat)
-          if (iostat == 0) then
+          IF (iostat == 0) THEN
              READ(98,*,end=998) nCFMIP
 998          CONTINUE
@@ -1426,6 +1425,6 @@
              IF(nCFMIP.GT.npCFMIP) THEN
                 print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
-                call abort_physic("physiq", "", 1)
-             else
+                CALL abort_physic("physiq", "", 1)
+             ELSE
                 print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
              ENDIF
@@ -1454,10 +1453,10 @@
                   tabijGCM, lonGCM, latGCM, iGCM, jGCM)
              !
-          else
+          ELSE
              ALLOCATE(tabijGCM(0))
              ALLOCATE(lonGCM(0), latGCM(0))
              ALLOCATE(iGCM(0), jGCM(0))
-          end if
-       else
+          ENDIF
+       ELSE
           ALLOCATE(tabijGCM(0))
           ALLOCATE(lonGCM(0), latGCM(0))
@@ -1490,8 +1489,8 @@
              zuthe(i)=0.
              zvthe(i)=0.
-             if(zstd(i).gt.10.)then
+             IF (zstd(i).gt.10.) THEN
                 zuthe(i)=(1.-zgam(i))*cos(zthe(i))
                 zvthe(i)=(1.-zgam(i))*sin(zthe(i))
-             endif
+             ENDIF
           ENDDO
        ENDIF
@@ -1534,5 +1533,5 @@
        ok_sync_omp=.false.
        CALL getin('ok_sync',ok_sync_omp)
-       call phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
+       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
             iGCM,jGCM,lonGCM,latGCM, &
             jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
@@ -1608,5 +1607,5 @@
           CALL VTb(VTphysiq)
 #endif
-       END IF
+       ENDIF
        !
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1614,10 +1613,10 @@
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-       call iniradia(klon,klev,paprs(1,1:klev+1))
+       CALL iniradia(klon,klev,paprs(1,1:klev+1))
 
        !$omp single
-       if (read_climoz >= 1) then
-          call open_climoz(ncid_climoz, press_climoz)
-       END IF
+       IF (read_climoz >= 1) THEN
+          CALL open_climoz(ncid_climoz, press_climoz)
+       ENDIF
        !$omp end single
        !
@@ -1700,5 +1699,5 @@
        CALL Rtime(debut)
 #endif
-    END IF
+    ENDIF
 
 
@@ -1760,9 +1759,9 @@
           ql_seri(i,k) = qx(i,k,iliq)
           !CR: ATTENTION, on rajoute la variable glace
-          if (nqo.eq.2) then
+          IF (nqo.eq.2) THEN
              qs_seri(i,k) = 0.
-          else if (nqo.eq.3) then
+          ELSE IF (nqo.eq.3) THEN
              qs_seri(i,k) = qx(i,k,isol)
-          endif
+          ENDIF
        ENDDO
     ENDDO
@@ -1808,5 +1807,5 @@
     ENDDO
     ! Initialize variables used for diagnostic purpose
-    if (flag_inhib_tend .ne. 0) call init_cmp_seri
+    IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
     !IM
     IF (ip_ebil_phy.ge.1) THEN 
@@ -1819,10 +1818,10 @@
        !     est egale a la variation de la physique au pas de temps precedent.
        !     Donc la somme de ces 2 variations devrait etre nulle.
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
 
     ! Diagnostiquer la tendance dynamique
@@ -1935,17 +1934,17 @@
       ELSE
         ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1   
-        if (ro3i == 361) ro3i = 360
-        if (read_climoz == 1) then
-           call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
+        IF (ro3i == 361) ro3i = 360
+        IF (read_climoz == 1) THEN
+           CALL regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
                 press_in_edg=press_climoz, paprs=paprs, v3=wo)
-        else
+        ELSE
            ! read_climoz == 2
-           call regr_pr_av(ncid_climoz, (/"tro3         ", &
+           CALL regr_pr_av(ncid_climoz, (/"tro3         ", &
                 "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
                 paprs=paprs, v3=wo)
-        end if
+        ENDIF
         ! Convert from mole fraction of ozone to column density of ozone in a
         ! cell, in kDU:
-        forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
+        FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
              * zmasse / dobson_u / 1e3
         ! (By regridding ozone values for LMDZ only once every 360th of
@@ -1971,5 +1970,5 @@
           !>jyg
 
-          if (iflag_ice_thermo.eq.0) then   
+          IF (iflag_ice_thermo.eq.0) THEN   
              !pas necessaire a priori
 
@@ -1984,5 +1983,5 @@
              d_q_eva(i,k) = zb
 
-          else
+          ELSE
 
              !CR: on r\'e-\'evapore eau liquide et glace
@@ -2002,5 +2001,5 @@
              d_t_eva(i,k) = za
              d_q_eva(i,k) = zb
-          endif
+          ENDIF
 
        ENDDO
@@ -2012,5 +2011,5 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
@@ -2018,5 +2017,5 @@
             , fs_bound, fq_bound )
        !
-    END IF
+    ENDIF
 
     !
@@ -2027,23 +2026,23 @@
 
     ! !!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
+    CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
     day_since_equinox = (jD_cur + jH_cur) - jD_eq
     !
     !   choix entre calcul de la longitude solaire vraie ou valeur fixee a 
     !   solarlong0
-    if (solarlong0<-999.) then
-       if (new_orbit) then
+    IF (solarlong0<-999.) THEN
+       IF (new_orbit) THEN
           ! calcul selon la routine utilisee pour les planetes
-          call solarlong(day_since_equinox, zlongi, dist)
-       else
+          CALL solarlong(day_since_equinox, zlongi, dist)
+       ELSE
           ! calcul selon la routine utilisee pour l'AR4
           CALL orbite(REAL(days_elapsed+1),zlongi,dist)
-       endif
-    else
+       ENDIF
+    ELSE
        zlongi=solarlong0  ! longitude solaire vraie
        dist=1.            ! distance au soleil / moyenne 
-    endif
-    if(prt_level.ge.1)                                                &
-         write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
+    ENDIF
+
+    IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
 
 
@@ -2055,6 +2054,6 @@
     ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
     ! non nul aux poles.
-    IF (abs(solarlong0-1000.)<1.e-4) then
-       call zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
+    IF (abs(solarlong0-1000.)<1.e-4) THEN
+       CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
             latitude_deg,longitude_deg,rmu0,fract)
        JrNt = 1.0
@@ -2109,10 +2108,10 @@
     ENDIF
 
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -2140,5 +2139,5 @@
 
 
-    if (iflag_pbl/=0) then
+    IF (iflag_pbl/=0) THEN
 
        !jyg+nrlmd<
@@ -2249,11 +2248,10 @@
        !--------------------------------------------------------------------
 
-       if (mydebug) then
-          call writefield_phy('u_seri',u_seri,nbp_lev)
-          call writefield_phy('v_seri',v_seri,nbp_lev)
-          call writefield_phy('t_seri',t_seri,nbp_lev)
-          call writefield_phy('q_seri',q_seri,nbp_lev)
-       endif
-
+       IF (mydebug) THEN
+          CALL writefield_phy('u_seri',u_seri,nbp_lev)
+          CALL writefield_phy('v_seri',v_seri,nbp_lev)
+          CALL writefield_phy('t_seri',t_seri,nbp_lev)
+          CALL writefield_phy('q_seri',q_seri,nbp_lev)
+       ENDIF
 
        !albedo SB >>>
@@ -2262,11 +2260,11 @@
        falb1=0.
        falb2=0.
-       select case(nsw)
-       case(2)
+       SELECT CASE(nsw)
+       CASE(2)
           albsol1=albsol_dir(:,1)
           albsol2=albsol_dir(:,2)
           falb1=falb_dir(:,1,:)
           falb2=falb_dir(:,2,:)
-       case(4)
+       CASE(4)
           albsol1=albsol_dir(:,1)
           albsol2=albsol_dir(:,2)*SFRWL(2)+albsol_dir(:,3)*SFRWL(3) &
@@ -2277,5 +2275,5 @@
                +falb_dir(:,4,:)*SFRWL(4)
           falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
-       case(6)
+       CASE(6)
           albsol1=albsol_dir(:,1)*SFRWL(1)+albsol_dir(:,2)*SFRWL(2) &
                +albsol_dir(:,3)*SFRWL(3)
@@ -2290,5 +2288,5 @@
                +falb_dir(:,6,:)*SFRWL(6)
           falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
-       end select
+       END SELECt
        !albedo SB <<<
 
@@ -2303,10 +2301,10 @@
                , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
                , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-          call diagphy(cell_area,ztit,ip_ebil_phy &
+          CALL diagphy(cell_area,ztit,ip_ebil_phy &
                , zero_v, zero_v, zero_v, zero_v, sens &
                , evap  , zero_v, zero_v, ztsol &
                , d_h_vcol, d_qt, d_ec &
                , fs_bound, fq_bound )
-       END IF
+       ENDIF
 
     ENDIF
@@ -2335,8 +2333,8 @@
     ENDDO
 
-    if (prt_level.ge.1) then
+    IF (prt_level.ge.1) THEN
        write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
        write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
-    endif
+    ENDIF
     !
     ! Appeler la convection (au choix)
@@ -2372,7 +2370,8 @@
        DO i = 1, klon
           omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
-       END DO
-    END DO
-    if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
+       ENDDO
+    ENDDO
+
+    IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
          omega(igout, :)
 
@@ -2408,23 +2407,19 @@
        !ajout pour la parametrisation des poches froides: calcul de
        !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
-       if (iflag_wake>=1) then
-         do k=1,klev
-            do i=1,klon
-                t_w(i,k) = t_seri(i,k) &
-                     +(1-wake_s(i))*wake_deltat(i,k)
-                q_w(i,k) = q_seri(i,k) &
-                     +(1-wake_s(i))*wake_deltaq(i,k)
-                t_x(i,k) = t_seri(i,k) &
-                     -wake_s(i)*wake_deltat(i,k)
-                q_x(i,k) = q_seri(i,k) &
-                     -wake_s(i)*wake_deltaq(i,k)
-            enddo
-         enddo
-       else
-                t_w(:,:) = t_seri(:,:)
+       IF (iflag_wake>=1) THEN
+         DO k=1,klev
+            DO i=1,klon
+                t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
+                q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
+                t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
+                q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
+            ENDDO
+         ENDDO
+       ELSE
+               t_w(:,:) = t_seri(:,:)
                 q_w(:,:) = q_seri(:,:)
                 t_x(:,:) = t_seri(:,:)
                 q_x(:,:) = q_seri(:,:)
-       endif
+       ENDIF
        !
        !jyg<
@@ -2495,5 +2490,5 @@
           ELSE
              nbtr_tmp=nbtr
-          END IF
+          ENDIF
           !jyg   iflag_con est dans clesphys
           !c          CALL concvl (iflag_con,iflag_clos,
@@ -2528,7 +2523,7 @@
           pmfu(:,:)=upwd(:,:)+dnwd(:,:)
 
-          do i = 1, klon
-             if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
-          enddo
+          DO i = 1, klon
+             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
+          ENDDO
           !
           !jyg<
@@ -2582,8 +2577,8 @@
        clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
        IF (iflag_cld_cv == 0) THEN
-          call clouds_gno &
+          CALL clouds_gno &
                (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
        ELSE
-          call clouds_bigauss &
+          CALL clouds_bigauss &
                (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
        ENDIF
@@ -2605,10 +2600,10 @@
           ema_pct(i)  = paprs(i,itop_con(i)+1)
 
-          if (itop_con(i).gt.klev-3) then
-             if(prt_level >= 9) then
+          IF (itop_con(i).gt.klev-3) THEN
+             IF (prt_level >= 9) THEN
                 write(lunout,*)'La convection monte trop haut '
                 write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
-             endif
-          endif
+             ENDIF
+          ENDIF
        ENDDO
     ELSE IF (iflag_con.eq.0) THEN
@@ -2626,5 +2621,5 @@
     ELSE
        WRITE(lunout,*) "iflag_con non-prevu", iflag_con
-       call abort_physic("physiq", "", 1)
+       CALL abort_physic("physiq", "", 1)
     ENDIF
 
@@ -2637,10 +2632,10 @@
     !-------------------------------------------------------------------------
 
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     !IM
@@ -2650,10 +2645,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, rain_con, snow_con, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
     !
     IF (check) THEN
@@ -2705,5 +2700,5 @@
     ! froides
     !
-    if (iflag_wake>=1) then
+    IF (iflag_wake>=1) THEN
        DO k=1,klev
           DO i=1,klon
@@ -2744,5 +2739,5 @@
        !
        !calcul caracteristiques de la poche froide
-       call calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
+       CALL calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
             t_seri, q_seri, omega,  &
             dt_dwn, dq_dwn, M_dwn, M_up,  &
@@ -2775,5 +2770,5 @@
        ENDIF   ! (iflag_wake_tend .GT. 0.)
 
-    endif  ! (iflag_wake>=1)
+    ENDIF  ! (iflag_wake>=1)
     !
     !===================================================================
@@ -2784,10 +2779,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
 
     !      print*,'apres callwake iflag_cld_th=', iflag_cld_th
@@ -2797,5 +2792,5 @@
     !===================================================================
     !
-    call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
+    CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
          ,seuil_inversion,weak_inversion,dthmin)
 
@@ -2814,18 +2809,18 @@
     !      detr_therm(:,:)=0.
     !
-    IF(prt_level>9)WRITE(lunout,*) &
+    IF (prt_level>9) WRITE(lunout,*) &
          'AVANT LA CONVECTION SECHE , iflag_thermals=' &
          ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
-    if(iflag_thermals<0) then
+    IF (iflag_thermals<0) THEN
        !  Rien
        !  ====
-       IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
-
-
-    else
+       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
+
+
+    ELSE
 
        !  Thermiques
        !  ==========
-       IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
+       IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
             ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
 
@@ -2842,5 +2837,5 @@
        !cc fin nrlmd le 10/04/2012
 
-       if (iflag_thermals>=1) then
+       IF (iflag_thermals>=1) THEN
           !jyg<
           IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
@@ -2866,5 +2861,5 @@
           ENDIF
           !>jyg
-          call calltherm(pdtphys &
+          CALL calltherm(pdtphys &
                ,pplay,paprs,pphi,weak_inversion &
                         ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
@@ -2932,12 +2927,11 @@
           ! -------------------------------------------------------------------
 
-          do i=1,klon
+          DO i=1,klon
              !           zmax_th(i)=pphi(i,lmax_th(i))/rg
              !CR:04/05/12:correction calcul zmax
              zmax_th(i)=zmax0(i) 
-          enddo
-
-       endif
-
+          ENDDO
+
+       ENDIF
 
        !  Ajustement sec
@@ -2948,12 +2942,12 @@
        ! Dans le cas contraire, on demarre au niveau 1.
 
-       if (iflag_thermals>=13.or.iflag_thermals<=0) then
-
-          if(iflag_thermals.eq.0) then
-             IF(prt_level>9)WRITE(lunout,*)'ajsec'
+       IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
+
+          IF (iflag_thermals.eq.0) THEN
+             IF (prt_level>9) WRITE(lunout,*)'ajsec'
              limbas(:)=1
-          else
+          ELSE
              limbas(:)=lmax_th(:)
-          endif
+          ENDIF
 
           ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
@@ -2963,13 +2957,13 @@
           ! non nulles numeriquement pour des mailles non concernees.
 
-          if (iflag_thermals==0) then
+          IF (iflag_thermals==0) THEN
              ! Calling adjustment alone (but not the thermal plume model)
              CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
                   , d_t_ajsb, d_q_ajsb)
-          else if (iflag_thermals>0) then
+          ELSE IF (iflag_thermals>0) THEN
              ! Calling adjustment above the top of thermal plumes
              CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
                   , d_t_ajsb, d_q_ajsb)
-          endif
+          ENDIF
 
           !--------------------------------------------------------------------
@@ -2982,7 +2976,7 @@
           !---------------------------------------------------------------------
 
-       endif
-
-    endif
+       ENDIF
+
+    ENDIF
     !
     !===================================================================
@@ -2993,10 +2987,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
 
 
@@ -3083,17 +3077,17 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, rain_lsc, snow_lsc, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
-
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    ENDIF
+
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     !
@@ -3110,19 +3104,19 @@
        !     print*,'avant calcul de la pseudo precip '
        !     print*,'iflag_cld_th',iflag_cld_th
-       if (iflag_cld_th.eq.-1) then
+       IF (iflag_cld_th.eq.-1) THEN
           rain_tiedtke=rain_con
-       else
+       ELSE
           !       print*,'calcul de la pseudo precip '
           rain_tiedtke=0.
           !         print*,'calcul de la pseudo precip 0'
-          do k=1,klev
-             do i=1,klon
-                if (d_q_con(i,k).lt.0.) then
+          DO k=1,klev
+             DO i=1,klon
+                IF (d_q_con(i,k).lt.0.) THEN
                    rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
                         *(paprs(i,k)-paprs(i,k+1))/rg
-                endif
-             enddo
-          enddo
-       endif
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDIF
        !
        !     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
@@ -3148,49 +3142,48 @@
        !  facttemps
        facteur = pdtphys *facttemps
-       do k=1,klev
-          do i=1,klon
+       DO k=1,klev
+          DO i=1,klon
              rnebcon(i,k)=rnebcon(i,k)*facteur
-             if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k)) &
-                  then
+             IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
                 rnebcon(i,k)=rnebcon0(i,k)
                 clwcon(i,k)=clwcon0(i,k)
-             endif
-          enddo
-       enddo
+             ENDIF
+          ENDDO
+       ENDDO
 
        !   On prend la somme des fractions nuageuses et des contenus en eau
 
-       if (iflag_cld_th>=5) then
-
-          do k=1,klev
+       IF (iflag_cld_th>=5) THEN
+
+          DO k=1,klev
              ptconvth(:,k)=fm_therm(:,k+1)>0.
-          enddo
-
-          if (iflag_coupl==4) then
+          ENDDO
+
+          IF (iflag_coupl==4) THEN
 
              ! Dans le cas iflag_coupl==4, on prend la somme des convertures
              ! convectives et lsc dans la partie des thermiques
              ! Le controle par iflag_coupl est peut etre provisoire.
-             do k=1,klev
-                do i=1,klon
-                   if (ptconv(i,k).and.ptconvth(i,k)) then
+             DO k=1,klev
+                DO i=1,klon
+                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
                       cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
                       cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
-                   else if (ptconv(i,k)) then
+                   ELSE IF (ptconv(i,k)) THEN
                       cldfra(i,k)=rnebcon(i,k)
                       cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
-                   endif
-                enddo
-             enddo
-
-          else if (iflag_coupl==5) then
-             do k=1,klev
-                do i=1,klon
+                   ENDIF
+                ENDDO
+             ENDDO
+
+          ELSE IF (iflag_coupl==5) THEN
+             DO k=1,klev
+                DO i=1,klon
                    cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
                    cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
-                enddo
-             enddo
-
-          else
+                ENDDO
+             ENDDO
+
+          ELSE
 
              ! Si on est sur un point touche par la convection
@@ -3202,21 +3195,21 @@
              ! definition des points sur lesquels ls thermiques sont actifs
 
-             do k=1,klev
-                do i=1,klon
-                   if (ptconv(i,k).and. .not. ptconvth(i,k)) then
+             DO k=1,klev
+                DO i=1,klon
+                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
                       cldfra(i,k)=rnebcon(i,k)
                       cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
-                   endif
-                enddo
-             enddo
-
-          endif
-
-       else
+                   ENDIF
+                ENDDO
+             ENDDO
+
+          ENDIF
+
+       ELSE
 
           ! Ancienne version
           cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
           cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
-       endif
+       ENDIF
 
     ENDIF
@@ -3256,10 +3249,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
     !
     ! Calculer l'humidite relative pour diagnostique
@@ -3319,5 +3312,5 @@
        calday = REAL(days_elapsed + 1) + jH_cur
 
-       call chemtime(itap+itau_phy-1, date0, dtime, itap)
+       CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
        IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
           CALL AEROSOL_METEO_CALC( &
@@ -3325,5 +3318,5 @@
                prfl,psfl,pctsrf,cell_area, &
                latitude_deg,longitude_deg,u10m,v10m)
-       END IF
+       ENDIF
 
        zxsnow_dummy(:) = 0.0
@@ -3368,5 +3361,5 @@
        CALL VTb(VTphysiq)
 #endif 
-    END IF !type_trac = inca
+    ENDIF !type_trac = inca
 
 
@@ -3393,5 +3386,5 @@
              IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
                 abort_message='config_inca=aero et rrtm=1 impossible'
-                call abort_physic(modname,abort_message,1)
+                CALL abort_physic(modname,abort_message,1)
              ELSE
                 !
@@ -3428,5 +3421,5 @@
                    abort_message='Only NSW=2 or 6 are possible with ' &
                         // 'aerosols and iflag_rrtm=1'
-                   call abort_physic(modname,abort_message,1)
+                   CALL abort_physic(modname,abort_message,1)
                 ENDIF
 
@@ -3438,5 +3431,5 @@
                 abort_message='You should compile with -rrtm if running ' &
                      // 'with iflag_rrtm=1'
-                call abort_physic(modname,abort_message,1)
+                CALL abort_physic(modname,abort_message,1)
 #endif
                 !
@@ -3497,4 +3490,6 @@
 #ifdef CPP_RRTM
 #ifdef CPP_StratAer
+       !--compute stratospheric mask
+       CALL stratosphere_mask(t_seri, pplay, latitude_deg)
        !--interactive strat aerosols
        CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
@@ -3510,5 +3505,5 @@
           mass_solu_aero(:,:)    = ccm(:,:,1) 
           mass_solu_aero_pi(:,:) = ccm(:,:,2) 
-       END IF
+       ENDIF
 
        IF (ok_newmicro) then
@@ -3619,10 +3614,10 @@
        ENDIF
 
-       if (mydebug) then
-          call writefield_phy('u_seri',u_seri,nbp_lev)
-          call writefield_phy('v_seri',v_seri,nbp_lev)
-          call writefield_phy('t_seri',t_seri,nbp_lev)
-          call writefield_phy('q_seri',q_seri,nbp_lev)
-       endif
+       IF (mydebug) THEN
+          CALL writefield_phy('u_seri',u_seri,nbp_lev)
+          CALL writefield_phy('v_seri',v_seri,nbp_lev)
+          CALL writefield_phy('t_seri',t_seri,nbp_lev)
+          CALL writefield_phy('q_seri',q_seri,nbp_lev)
+       ENDIF
 
        !
@@ -3632,5 +3627,5 @@
        IF (iflag_radia .ge. 2) THEN
           zsav_tsol (:) = zxtsol(:)
-          call perturb_radlwsw(zxtsol,iflag_radia)
+          CALL perturb_radlwsw(zxtsol,iflag_radia)
        ENDIF
 
@@ -3721,8 +3716,8 @@
           !IM Par defaut on a les taux perturbes egaux aux taux actuels
           !
-          if (ok_4xCO2atm) then
-             if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
-                  RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
-                  RCFC12_per.NE.RCFC12_act) THEN
+          IF (ok_4xCO2atm) THEN
+             IF (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.     &
+                 RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
+                 RCFC12_per.NE.RCFC12_act) THEN
                 !
                 RCO2 = RCO2_per
@@ -3792,5 +3787,5 @@
           PRINT *,'>>>>           heat et cool mis a zero '
           PRINT *,'--------------------------------------------------'
-       END IF
+       ENDIF
        heat=0.
        cool=0.
@@ -3804,5 +3799,5 @@
        lwdn=0.
        lwdn0=0.
-    END IF
+    ENDIF
 
     !
@@ -3813,7 +3808,7 @@
     radsol=solsw*swradcorr+sollw
 
-    if (ok_4xCO2atm) then
+    IF (ok_4xCO2atm) THEN
        radsolp=solswp*swradcorr+sollwp
-    endif
+    ENDIF
 
     !
@@ -3833,10 +3828,10 @@
 
     !
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     !IM
@@ -3846,10 +3841,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , topsw, toplw, solsw, sollw, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
     !
     !
@@ -3919,10 +3914,10 @@
     ENDIF ! fin de test sur ok_orodr
     !
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     IF (ok_orolf) THEN
@@ -4003,6 +3998,6 @@
     ENDIF
 
-    if (ok_gwd_rando) then
-       call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
+    IF (ok_gwd_rando) THEN
+       CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
             rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
             du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
@@ -4017,14 +4012,14 @@
                * (paprs(:, k)-paprs(:, k+1))/rg
        ENDDO
-    end if
+    ENDIF
 
     ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
 
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     DO i = 1, klon
@@ -4059,10 +4054,10 @@
             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , zero_v, zero_v, zero_v, zero_v, zero_v &
             , zero_v, zero_v, zero_v, ztsol &
             , d_h_vcol, d_qt, d_ec &
             , fs_bound, fq_bound )
-    END IF
+    ENDIF
 
     !DC Calcul de la tendance due au methane
@@ -4072,5 +4067,5 @@
        CALL add_phys_tend(du0, dv0, dt0, d_q_ch4*dtime, dql0, dqi0, paprs, &
             'q_ch4', abortphy,flag_inhib_tend)
-    END IF
+    ENDIF
     !
     !
@@ -4091,5 +4086,5 @@
           !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
           !     s        ref_liq,ref_ice
-          call phys_cosp(itap,dtime,freq_cosp, &
+          CALL phys_cosp(itap,dtime,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
                ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
@@ -4121,16 +4116,15 @@
 
   IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
-  write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', &
-     & ok_airs, freq_airs
-  call simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
-     & map_prop_hc,map_prop_hist,&
-     & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
-     & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
-     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
-     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
-     & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
-     & map_ntot,map_hc,map_hist,&
-     & map_Cb,map_ThCi,map_Anv,&
-     & alt_tropo )
+     write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
+     CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
+        & map_prop_hc,map_prop_hist,&
+        & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
+        & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
+        & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
+        & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
+        & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
+        & map_ntot,map_hc,map_hist,&
+        & map_Cb,map_ThCi,map_Anv,&
+        & alt_tropo )
   ENDIF
 
@@ -4151,5 +4145,5 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-    END IF
+    ENDIF
 
 #ifdef CPP_Dust
@@ -4171,5 +4165,5 @@
 #else
 
-    call phytrac ( &
+    CALL phytrac ( &
          itap,     days_elapsed+1,    jH_cur,   debut, &
          lafin,    dtime,     u, v,     t, &
@@ -4201,5 +4195,5 @@
        IF (prt_level.ge.9) &
             print*,'Attention on met a 0 les thermiques pour phystoke'
-       call phystokenc ( &
+       CALL phystokenc ( &
             nlon,klev,pdtphys,longitude_deg,latitude_deg, &
             t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
@@ -4255,5 +4249,5 @@
        !     Donc la somme de ces 2 variations devrait etre nulle.
 
-       call diagphy(cell_area,ztit,ip_ebil_phy &
+       CALL diagphy(cell_area,ztit,ip_ebil_phy &
             , topsw, toplw, solsw, sollw, sens &
             , evap, rain_fall, snow_fall, ztsol &
@@ -4263,5 +4257,5 @@
        d_h_vcol_phy=d_h_vcol
        !
-    END IF
+    ENDIF
     !
     !=======================================================================
@@ -4326,5 +4320,5 @@
        CALL VTb(VTphysiq)
 #endif
-    END IF
+    ENDIF
 
 
@@ -4336,10 +4330,10 @@
     ENDIF
     !
-    if (mydebug) then
-       call writefield_phy('u_seri',u_seri,nbp_lev)
-       call writefield_phy('v_seri',v_seri,nbp_lev)
-       call writefield_phy('t_seri',t_seri,nbp_lev)
-       call writefield_phy('q_seri',q_seri,nbp_lev)
-    endif
+    IF (mydebug) THEN
+       CALL writefield_phy('u_seri',u_seri,nbp_lev)
+       CALL writefield_phy('v_seri',v_seri,nbp_lev)
+       CALL writefield_phy('t_seri',t_seri,nbp_lev)
+       CALL writefield_phy('q_seri',q_seri,nbp_lev)
+    ENDIF
 
     DO k = 1, klev
@@ -4351,7 +4345,7 @@
           d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
           !CR: on ajoute le contenu en glace
-          if (nqo.eq.3) then
+          IF (nqo.eq.3) THEN
              d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
-          endif
+          ENDIF
        ENDDO
     ENDDO
@@ -4420,5 +4414,5 @@
     !==========================================================================
 
-    if (prt_level.ge.1) then
+    IF (prt_level.ge.1) THEN
        write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
        write(lunout,*) &
@@ -4429,13 +4423,13 @@
             pctsrf(igout,is_sic)
        write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
                d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
                d_t_eva(igout,k)
-       enddo
+       ENDDO
        write(lunout,*) 'cool,heat'
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) cool(igout,k),heat(igout,k)
-       enddo
+       ENDDO
 
        !jyg<     (En attendant de statuer sur le sort de d_t_oli)
@@ -4446,19 +4440,17 @@
        !jyg!     enddo
        write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) d_t_vdf(igout,k), &
                d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
-       enddo
+       ENDDO
        !>jyg
 
        write(lunout,*) 'd_ps ',d_ps(igout)
        write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
-       do k=1,klev
+       DO k=1,klev
           write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
                d_qx(igout,k,1),d_qx(igout,k,2)
-       enddo
-    endif
-
-    !==========================================================================
+       ENDDO
+    ENDIF
 
     !============================================================
@@ -4501,27 +4493,26 @@
     !=============================================================
 
-    if (iflag_thermals>=1) then
+    IF (iflag_thermals>=1) THEN
        d_t_lscth=0.
        d_t_lscst=0.
        d_q_lscth=0.
        d_q_lscst=0.
-       do k=1,klev
-          do i=1,klon
-             if (ptconvth(i,k)) then
+       DO k=1,klev
+          DO i=1,klon
+             IF (ptconvth(i,k)) THEN
                 d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
                 d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
-             else
+             ELSE
                 d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
                 d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
-             endif
-          enddo
-       enddo
-
-       do i=1,klon
+             ENDIF
+          ENDDO
+       ENDDO
+
+       DO i=1,klon
           plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
           plul_th(i)=prfl(i,1)+psfl(i,1)
-       enddo
-    endif
-
+       ENDDO
+    ENDIF
 
     !On effectue les sorties:
@@ -4543,5 +4534,4 @@
 #endif
 
-
 #ifndef CPP_XIOS
     CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
@@ -4575,10 +4565,10 @@
        !         close(97)
        !$OMP MASTER
-       if (read_climoz >= 1) then
-          if (is_mpi_root) then
-             call nf95_close(ncid_climoz)
-          end if
-          deallocate(press_climoz) ! pointer
-       end if
+       IF (read_climoz >= 1) THEN
+          IF (is_mpi_root) THEN
+             CALL nf95_close(ncid_climoz)
+          ENDIF
+          DEALLOCATE(press_climoz) ! pointer
+       ENDIF
        !$OMP END MASTER
     ENDIF
