Index: LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 3389)
+++ LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 3395)
@@ -7,7 +7,7 @@
 
 CONTAINS
-!
-!******************************************************************************
-!
+  !
+  !******************************************************************************
+  !
   SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
        windsp, rmu0, fder, tsurf_in, &
@@ -22,16 +22,16 @@
        flux_u1, flux_v1)
 
-  use albedo, only: alboc, alboc_cd
-  USE dimphy, ONLY: klon, zmasq
-  USE surface_data, ONLY     : type_ocean
-  USE ocean_forced_mod, ONLY : ocean_forced_noice
-  USE ocean_slab_mod, ONLY   : ocean_slab_noice
-  USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
-  USE indice_sol_mod, ONLY : nbsrf, is_oce
-  USE limit_read_mod
-!
-! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
-! slab or couple). The calculations of albedo and rugosity for the ocean surface are 
-! done in here because they are identical for the different modes of ocean. 
+    use albedo, only: alboc, alboc_cd
+    USE dimphy, ONLY: klon, zmasq
+    USE surface_data, ONLY     : type_ocean
+    USE ocean_forced_mod, ONLY : ocean_forced_noice
+    USE ocean_slab_mod, ONLY   : ocean_slab_noice
+    USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
+    USE indice_sol_mod, ONLY : nbsrf, is_oce
+    USE limit_read_mod
+    !
+    ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
+    ! slab or couple). The calculations of albedo and rugosity for the ocean surface are 
+    ! done in here because they are identical for the different modes of ocean. 
 
 
@@ -41,6 +41,6 @@
     ! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0)
 
-! Input variables
-!******************************************************************************
+    ! Input variables
+    !******************************************************************************
     INTEGER, INTENT(IN)                      :: itime, jour, knon
     INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
@@ -66,6 +66,6 @@
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
 
-! In/Output variables
-!******************************************************************************
+    ! In/Output variables
+    !******************************************************************************
     REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
     REAL, DIMENSION(klon), INTENT(INOUT)     :: qsurf
@@ -73,13 +73,13 @@
     REAL, DIMENSION(klon), INTENT(inOUT):: z0h
 
-! Output variables
-!******************************************************************************
+    ! Output variables
+    !******************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)       :: z0m
-!albedo SB >>>
-!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
-!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
+    !albedo SB >>>
+    !    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
+    !    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
     REAL, DIMENSION(6), INTENT(IN)          :: SFRWL 
     REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
-!albedo SB <<<     
+    !albedo SB <<<     
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
@@ -88,6 +88,6 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
 
-! Local variables
-!******************************************************************************
+    ! Local variables
+    !******************************************************************************
     INTEGER               :: i, k
     REAL                  :: tmp
@@ -98,28 +98,28 @@
     CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
 
-! End definition
-!******************************************************************************
-
-
-!******************************************************************************
-! Calculate total net radiance at surface
-!
-!******************************************************************************
+    ! End definition
+    !******************************************************************************
+
+
+    !******************************************************************************
+    ! Calculate total net radiance at surface
+    !
+    !******************************************************************************
     radsol(1:klon) = 0.0 ! initialisation a priori inutile
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
 
-!******************************************************************************
-! Cdragq computed from cdrag
-! The difference comes only from a factor (f_z0qh_oce) on z0, so that
-! it can be computed inside surf_ocean
-! More complicated appraches may require the propagation through
-! pbl_surface of an independant cdragq variable.
-!******************************************************************************
+    !******************************************************************************
+    ! Cdragq computed from cdrag
+    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
+    ! it can be computed inside surf_ocean
+    ! More complicated appraches may require the propagation through
+    ! pbl_surface of an independant cdragq variable.
+    !******************************************************************************
 
     IF ( f_z0qh_oce .ne. 1.) THEN
-! Si on suit les formulations par exemple de Tessel, on 
-! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
+       ! Si on suit les formulations par exemple de Tessel, on 
+       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
        cdragq(1:knon)=cdragh(1:knon)*                                      &
-       log(z1lay(1:knon)/z0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*z0h(1:knon)))
+            log(z1lay(1:knon)/z0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*z0h(1:knon)))
     ELSE
        cdragq(1:knon)=cdragh(1:knon)
@@ -127,7 +127,7 @@
 
 
-!******************************************************************************
-! Switch according to type of ocean (couple, slab or forced)
-!******************************************************************************
+    !******************************************************************************
+    ! Switch according to type of ocean (couple, slab or forced)
+    !******************************************************************************
     SELECT CASE(type_ocean)
     CASE('couple')
@@ -154,5 +154,5 @@
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
             tsurf_new, dflux_s, dflux_l, lmt_bils)
-       
+
     CASE('force')
        CALL ocean_forced_noice( &
@@ -168,108 +168,108 @@
     END SELECT
 
-!******************************************************************************
-! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
-!******************************************************************************
+    !******************************************************************************
+    ! fcodron: compute lmt_bils  forced case (same as wfbils_oce / 1.-contfracatm)
+    !******************************************************************************
     IF (type_ocean.NE.'slab') THEN
-        lmt_bils(1:klon)=0.
-        DO i=1,knon
-           lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
-           *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i))) 
-        END DO
+       lmt_bils(1:klon)=0.
+       DO i=1,knon
+          lmt_bils(knindex(i))=(swnet(i)+lwnet(i)+fluxsens(i)+fluxlat(i)) &
+               *pctsrf(knindex(i),is_oce)/(1.-zmasq(knindex(i))) 
+       END DO
     END IF
 
-!******************************************************************************
-! Calculate ocean surface albedo
-!******************************************************************************
-!albedo SB >>>
-IF (iflag_albedo==0) THEN
-!--old parametrizations of ocean surface albedo
-!
-    IF (iflag_cycle_diurne.GE.1) THEN
-!
-       CALL alboc_cd(rmu0,alb_eau)
-!
-!--ad-hoc correction for model radiative balance tuning
-!--now outside alboc_cd routine
-       alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
-       alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.0),1.0)
-!
+    !******************************************************************************
+    ! Calculate ocean surface albedo
+    !******************************************************************************
+    !albedo SB >>>
+    IF (iflag_albedo==0) THEN
+       !--old parametrizations of ocean surface albedo
+       !
+       IF (iflag_cycle_diurne.GE.1) THEN
+          !
+          CALL alboc_cd(rmu0,alb_eau)
+          !
+          !--ad-hoc correction for model radiative balance tuning
+          !--now outside alboc_cd routine
+          alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
+          alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.0),1.0)
+          !
+       ELSE
+          !
+          CALL alboc(REAL(jour),rlat,alb_eau)
+          !--ad-hoc correction for model radiative balance tuning
+          !--now outside alboc routine
+          alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
+          alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.04),0.60)
+          !
+       ENDIF
+       !
+       DO i =1, knon
+          DO  k=1,nsw
+             alb_dir_new(i,k) = alb_eau(knindex(i))
+          ENDDO
+       ENDDO
+       !IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions 
+       !albedo for diffuse radiation is taken the same as for direct radiation
+       alb_dif_new(1:knon,:)=alb_dir_new(1:knon,:)
+       !IM 09122015 end
+       !
+    ELSE IF (iflag_albedo==1) THEN
+       !--new parametrization of ocean surface albedo by Sunghye Baek
+       !--albedo for direct and diffuse radiation are different
+       !
+       CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
+       !
+       !--ad-hoc correction for model radiative balance tuning
+       alb_dir_new(1:knon,:) = fmagic*alb_dir_new(1:knon,:) + pmagic
+       alb_dif_new(1:knon,:) = fmagic*alb_dif_new(1:knon,:) + pmagic
+       alb_dir_new(1:knon,:)=MIN(MAX(alb_dir_new(1:knon,:),0.0),1.0)
+       alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
+       !
+    ELSE IF (iflag_albedo==2) THEN
+       ! F. Codron albedo read from limit.nc
+       CALL limit_read_rug_alb(itime, dtime, jour,&
+            knon, knindex, z0_lim, alb_eau)
+       DO i =1, knon
+          DO  k=1,nsw
+             alb_dir_new(i,k) = alb_eau(i)
+          ENDDO
+       ENDDO
+       alb_dif_new=alb_dir_new
+    ENDIF
+    !albedo SB <<<
+
+    !******************************************************************************
+    ! Calculate the rugosity
+    !******************************************************************************
+    IF (iflag_z0_oce==0) THEN
+       DO i = 1, knon
+          tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
+          z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
+               +  0.11*14e-6 / SQRT(cdragm(i) * tmp)
+          z0m(i) = MAX(1.5e-05,z0m(i))
+       ENDDO
+       z0h(1:knon)=z0m(1:knon) ! En attendant mieux
+
+    ELSE IF (iflag_z0_oce==1) THEN
+       DO i = 1, knon
+          tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
+          z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
+               + 0.11*14e-6 / SQRT(cdragm(i) * tmp)
+          z0m(i) = MAX(1.5e-05,z0m(i))
+          z0h(i)=0.4*14e-6 / SQRT(cdragm(i) * tmp)
+       ENDDO
+    ELSE IF (iflag_z0_oce==-1) THEN
+       DO i = 1, knon
+          z0m(i) = z0min
+          z0h(i) = z0min
+       ENDDO
     ELSE
-!
-       CALL alboc(REAL(jour),rlat,alb_eau)
-!--ad-hoc correction for model radiative balance tuning
-!--now outside alboc routine
-       alb_eau(1:klon) = fmagic*alb_eau(1:klon) + pmagic
-       alb_eau(1:klon)=MIN(MAX(alb_eau(1:klon),0.04),0.60)
-!
+       CALL abort_physic(modname,'version non prevue',1)
     ENDIF
-!
-    DO i =1, knon
-      DO  k=1,nsw
-       alb_dir_new(i,k) = alb_eau(knindex(i))
-      ENDDO
-    ENDDO
-!IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions 
-!albedo for diffuse radiation is taken the same as for direct radiation
-     alb_dif_new(1:knon,:)=alb_dir_new(1:knon,:)
-!IM 09122015 end
-!
-ELSE IF (iflag_albedo==1) THEN
-!--new parametrization of ocean surface albedo by Sunghye Baek
-!--albedo for direct and diffuse radiation are different
-!
-    CALL ocean_albedo(knon,rmu0,knindex,windsp,SFRWL,alb_dir_new,alb_dif_new)
-!
-!--ad-hoc correction for model radiative balance tuning
-    alb_dir_new(1:knon,:) = fmagic*alb_dir_new(1:knon,:) + pmagic
-    alb_dif_new(1:knon,:) = fmagic*alb_dif_new(1:knon,:) + pmagic
-    alb_dir_new(1:knon,:)=MIN(MAX(alb_dir_new(1:knon,:),0.0),1.0)
-    alb_dif_new(1:knon,:)=MIN(MAX(alb_dif_new(1:knon,:),0.0),1.0)
-!
-! F. Codron albedo read from limit.nc
-ELSE IF (iflag_albedo==2) THEN
-    CALL limit_read_rug_alb(itime, dtime, jour,&
-         knon, knindex, z0_lim, alb_eau)
-    DO i =1, knon
-      DO  k=1,nsw
-       alb_dir_new(i,k) = alb_eau(i)
-      ENDDO
-    ENDDO
-    alb_dif_new=alb_dir_new
-ENDIF
-!albedo SB <<<
-
-!******************************************************************************
-! Calculate the rugosity
-!******************************************************************************
-IF (iflag_z0_oce==0) THEN
-    DO i = 1, knon
-       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
-       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
-            +  0.11*14e-6 / SQRT(cdragm(i) * tmp)
-       z0m(i) = MAX(1.5e-05,z0m(i))
-    ENDDO   
-    z0h(1:knon)=z0m(1:knon) ! En attendant mieux
-
-ELSE IF (iflag_z0_oce==1) THEN
-    DO i = 1, knon
-       tmp = MAX(cepdu2,gustiness(i)+u1(i)**2+v1(i)**2)
-       z0m(i) = 0.018*cdragm(i) * (gustiness(i)+u1(i)**2+v1(i)**2)/RG  &
-            + 0.11*14e-6 / SQRT(cdragm(i) * tmp)
-       z0m(i) = MAX(1.5e-05,z0m(i))
-       z0h(i)=0.4*14e-6 / SQRT(cdragm(i) * tmp)
-    ENDDO
-ELSE IF (iflag_z0_oce==-1) THEN
-    DO i = 1, knon
-       z0m(i) = z0min
-       z0h(i) = z0min
-    ENDDO
-ELSE
-       CALL abort_physic(modname,'version non prevue',1)
-ENDIF
-!
-!******************************************************************************
+    !
+    !******************************************************************************
   END SUBROUTINE surf_ocean
-!******************************************************************************
-!
+  !******************************************************************************
+  !
 END MODULE surf_ocean_mod
