Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cpl_oasis3.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cpl_oasis3.F90	(revision 11963)
@@ -114,5 +114,5 @@
       !------------------------------------------------------------------
       CALL oasis_init_comp ( ncomp_id, TRIM(cd_modname), nerror )
-      IF ( nerror /= OASIS_Ok ) &
+      IF( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init', 'Failure in oasis_init_comp')
 
@@ -122,5 +122,5 @@
 
       CALL oasis_get_localcomm ( kl_comm, nerror )
-      IF ( nerror /= OASIS_Ok ) &
+      IF( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
       !
@@ -149,5 +149,5 @@
 
       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define
-      IF ( ltmp_wapatch ) THEN
+      IF( ltmp_wapatch ) THEN
          nldi_save = nldi   ;   nlei_save = nlei
          nldj_save = nldj   ;   nlej_save = nlej
@@ -217,5 +217,5 @@
       !
       DO ji = 1, ksnd
-         IF ( ssnd(ji)%laction ) THEN
+         IF( ssnd(ji)%laction ) THEN
 
             IF( ssnd(ji)%nct > nmaxcat ) THEN
@@ -228,5 +228,5 @@
                DO jm = 1, kcplmodel
 
-                  IF ( ssnd(ji)%nct .GT. 1 ) THEN
+                  IF( ssnd(ji)%nct .GT. 1 ) THEN
                      WRITE(cli2,'(i2.2)') jc
                      zclname = TRIM(ssnd(ji)%clname)//'_cat'//cli2
@@ -234,5 +234,5 @@
                      zclname = ssnd(ji)%clname
                   ENDIF
-                  IF ( kcplmodel  > 1 ) THEN
+                  IF( kcplmodel  > 1 ) THEN
                      WRITE(cli2,'(i2.2)') jm
                      zclname = 'model'//cli2//'_'//TRIM(zclname)
@@ -241,10 +241,10 @@
                   IF( agrif_fixed() /= 0 ) THEN 
                      zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
-                  END IF
+                  ENDIF
 #endif
                   IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
                   CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
                      &                OASIS_Out          , ishape , OASIS_REAL, nerror )
-                  IF ( nerror /= OASIS_Ok ) THEN
+                  IF( nerror /= OASIS_Ok ) THEN
                      WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
                      CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
@@ -262,5 +262,5 @@
       !
       DO ji = 1, krcv
-         IF ( srcv(ji)%laction ) THEN 
+         IF( srcv(ji)%laction ) THEN 
             
             IF( srcv(ji)%nct > nmaxcat ) THEN
@@ -273,5 +273,5 @@
                DO jm = 1, kcplmodel
                   
-                  IF ( srcv(ji)%nct .GT. 1 ) THEN
+                  IF( srcv(ji)%nct .GT. 1 ) THEN
                      WRITE(cli2,'(i2.2)') jc
                      zclname = TRIM(srcv(ji)%clname)//'_cat'//cli2
@@ -279,5 +279,5 @@
                      zclname = srcv(ji)%clname
                   ENDIF
-                  IF ( kcplmodel  > 1 ) THEN
+                  IF( kcplmodel  > 1 ) THEN
                      WRITE(cli2,'(i2.2)') jm
                      zclname = 'model'//cli2//'_'//TRIM(zclname)
@@ -286,10 +286,10 @@
                   IF( agrif_fixed() /= 0 ) THEN 
                      zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
-                  END IF
+                  ENDIF
 #endif
                   IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
                   CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
                      &                OASIS_In           , ishape , OASIS_REAL, nerror )
-                  IF ( nerror /= OASIS_Ok ) THEN
+                  IF( nerror /= OASIS_Ok ) THEN
                      WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
                      CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
@@ -310,5 +310,5 @@
       IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')
       !
-      IF ( ltmp_wapatch ) THEN
+      IF( ltmp_wapatch ) THEN
          nldi = nldi_save   ;   nlei = nlei_save
          nldj = nldj_save   ;   nlej = nlej_save
@@ -332,5 +332,5 @@
       !!--------------------------------------------------------------------
       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define
-      IF ( ltmp_wapatch ) THEN
+      IF( ltmp_wapatch ) THEN
          nldi_save = nldi   ;   nlei_save = nlei
          nldj_save = nldj   ;   nlej_save = nlej
@@ -349,6 +349,6 @@
                CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo )
                
-               IF ( ln_ctl ) THEN        
-                  IF ( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   &
+               IF( ln_ctl ) THEN        
+                  IF( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   &
                      & kinfo == OASIS_SentOut  .OR. kinfo == OASIS_ToRestOut ) THEN
                      WRITE(numout,*) '****************'
@@ -368,5 +368,5 @@
          ENDDO
       ENDDO
-      IF ( ltmp_wapatch ) THEN
+      IF( ltmp_wapatch ) THEN
          nldi = nldi_save   ;   nlei = nlei_save
          nldj = nldj_save   ;   nlej = nlej_save
@@ -393,5 +393,5 @@
       !!--------------------------------------------------------------------
       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define
-      IF ( ltmp_wapatch ) THEN
+      IF( ltmp_wapatch ) THEN
          nldi_save = nldi   ;   nlei_save = nlei
          nldj_save = nldj   ;   nlej_save = nlej
@@ -403,5 +403,5 @@
       !
       DO jc = 1, srcv(kid)%nct
-         IF ( ltmp_wapatch ) THEN
+         IF( ltmp_wapatch ) THEN
             IF( nimpp           ==      1 ) nldi = 1
             IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi
@@ -420,7 +420,7 @@
                   &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
                
-               IF ( ln_ctl )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
+               IF( ln_ctl )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
                
-               IF ( llaction ) THEN
+               IF( llaction ) THEN
                   
                   kinfo = OASIS_Rcv
@@ -432,5 +432,5 @@
                   ENDIF
                   
-                  IF ( ln_ctl ) THEN        
+                  IF( ln_ctl ) THEN        
                      WRITE(numout,*) '****************'
                      WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
@@ -450,5 +450,5 @@
          ENDDO
 
-         IF ( ltmp_wapatch ) THEN
+         IF( ltmp_wapatch ) THEN
             nldi = nldi_save   ;   nlei = nlei_save
             nldj = nldj_save   ;   nlej = nlej_save
@@ -483,5 +483,5 @@
       !
       DO ji = 1, nsnd
-         IF (ssnd(ji)%laction ) THEN
+         IF(ssnd(ji)%laction ) THEN
             DO jm = 1, ncplmodel
                IF( ssnd(ji)%nid(1,jm) /= -1 ) THEN
@@ -495,5 +495,5 @@
       ENDDO
       DO ji = 1, nrcv
-         IF (srcv(ji)%laction ) THEN
+         IF(srcv(ji)%laction ) THEN
             DO jm = 1, ncplmodel
                IF( srcv(ji)%nid(1,jm) /= -1 ) THEN
@@ -529,5 +529,5 @@
       !
       DEALLOCATE( exfld )
-      IF (nstop == 0) THEN
+      IF(nstop == 0) THEN
          CALL oasis_terminate( nerror )         
       ELSE
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cyclone.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cyclone.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cyclone.F90	(revision 11963)
@@ -137,5 +137,5 @@
             zhemi = SIGN( 1. , zrlat )
             zinfl = 15.* rad                             ! clim inflow angle in Tropical Cyclones
-         IF ( vortex == 0 ) THEN
+         IF( vortex == 0 ) THEN
 
             ! Vortex Holland reconstruct wind at each lon-lat position
@@ -157,19 +157,19 @@
                      &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) )
 
-                 IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius
+                 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius
                   ! shape of the wind profile
                   zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb
                   zztmp =  zvmax * SQRT( zztmp * EXP(1. - zztmp) )    
 
-                  IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2
+                  IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2
                      zztmp = zztmp * ( (zrout2-zdist)*1.e-6 )
                   ENDIF
 
                   ! !!! KILL EQ WINDS
-                  ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN
+                  ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN
                   !    zztmp = 0.                              ! winds in other hemisphere
-                  !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S
-                  ! ENDIF
-                  ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN
+                  !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S
+                  ! ENDIF
+                  ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN
                   !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 
                   !    !linear to zero between 10 and 5
@@ -177,5 +177,5 @@
                   ! !!! / KILL EQ
 
-                  IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude
+                  IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude
 
                   zwnd_t =   COS( zinfl ) * zztmp    
@@ -196,5 +196,5 @@
             END DO
          
-         ELSE IF ( vortex == 1 ) THEN
+         ELSE IF( vortex == 1 ) THEN
 
             ! Vortex Willoughby reconstruct wind at each lon-lat position
@@ -206,5 +206,5 @@
             zn   =   2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) )            
             zA   =   0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) )  
-            IF (zA < 0) THEN 
+            IF(zA < 0) THEN 
                zA=0
             ENDIF           
@@ -218,8 +218,8 @@
                      &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) )
 
-                 IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius
+                 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius
                
                   ! shape of the wind profile                     
-                  IF (zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind
+                  IF(zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind
                      zztmp  = zvmax * (zdist/zrmw)**zn
                   ELSE 
@@ -227,14 +227,14 @@
                   ENDIF
 
-                  IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2
+                  IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2
                      zztmp = zztmp * ( (zrout2-zdist)*1.e-6 )
                   ENDIF
 
                   ! !!! KILL EQ WINDS
-                  ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN
+                  ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN
                   !    zztmp = 0.                              ! winds in other hemisphere
-                  !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S
-                  ! ENDIF
-                  ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN
+                  !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S
+                  ! ENDIF
+                  ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN
                   !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 
                   !    !linear to zero between 10 and 5
@@ -242,5 +242,5 @@
                   ! !!! / KILL EQ
 
-                  IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude
+                  IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude
 
                   zwnd_t =   COS( zinfl ) * zztmp    
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/fldread.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/fldread.F90	(revision 11963)
@@ -167,5 +167,5 @@
       IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
 
-      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
+      IF( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
       ELSE                                      ;   it_offset = 0
       ENDIF
@@ -389,5 +389,5 @@
          ENDIF
          !
-         IF ( sdjf%cltype(1:4) == 'week' ) THEN
+         IF( sdjf%cltype(1:4) == 'week' ) THEN
             isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week
             llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month
@@ -464,5 +464,5 @@
       ENDIF
       !
-      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
+      IF( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
       ELSE                                      ;   it_offset = 0
       ENDIF
@@ -656,5 +656,5 @@
             ENDIF
          CASE DEFAULT
-            IF (lk_c1d .AND. lmoor ) THEN
+            IF(lk_c1d .AND. lmoor ) THEN
                IF( sdjf%ln_tint ) THEN
                   CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
@@ -1071,5 +1071,5 @@
          imonth = kmonth
          iday = kday
-         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
+         IF( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
             isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )  
             llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
@@ -1080,5 +1080,5 @@
          ENDIF
       ELSE                                                  ! use current day values
-         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
+         IF( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
             isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
             llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
@@ -1318,5 +1318,5 @@
 
       !! get dimensions
-      IF ( SIZE(sd%fnow, 3) > 1 ) THEN
+      IF( SIZE(sd%fnow, 3) > 1 ) THEN
          ALLOCATE( ddims(4) )
       ELSE
@@ -1331,5 +1331,5 @@
 
       CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
-      IF ( inum > 0 ) THEN
+      IF( inum > 0 ) THEN
 
          !! determine whether we have an east-west cyclic grid
@@ -1663,5 +1663,5 @@
       END DO
 
-      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
+      IF(ref_wgts(kw)%numwgt .EQ. 16) THEN
 
         !! fix up halo points that we couldnt read from file
@@ -1746,5 +1746,5 @@
          END DO
          !
-      END IF
+      ENDIF
       !
    END SUBROUTINE fld_interp
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcapr.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcapr.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcapr.F90	(revision 11963)
@@ -103,5 +103,5 @@
       !
       !                                            !* control check
-      IF ( ln_apr_obc  ) THEN
+      IF( ln_apr_obc  ) THEN
          IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data'
       ENDIF
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90	(revision 11963)
@@ -259,16 +259,16 @@
       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
       !
-      IF ( ln_wave ) THEN
+      IF( ln_wave ) THEN
          !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
+         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
             CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
             !drag coefficient read from wave model definable only with mfs bulk formulae and core
-         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN
+         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN
             CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
-         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
+         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN
             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
          ENDIF
       ELSE
-         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
+         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
             &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
             &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
@@ -481,8 +481,8 @@
       zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
 
-      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
          zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
-      END IF
+      ENDIF
 
       !LB:
@@ -492,8 +492,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
+         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
+         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
       END SELECT
@@ -506,5 +506,5 @@
 
 
-      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
 
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
@@ -543,5 +543,5 @@
          tsk(:,:) = zst(:,:)
 
-      ELSE !IF ( ln_skin_cs .OR. ln_skin_wl )
+      ELSE !IF( ln_skin_cs .OR. ln_skin_wl )
 
          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
@@ -567,5 +567,5 @@
          END SELECT
 
-      END IF ! IF ( ln_skin_cs .OR. ln_skin_wl )
+      ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl )
 
       !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
@@ -576,5 +576,5 @@
          t_zu(:,:) = ztpot(:,:)
          q_zu(:,:) = zqair(:,:)
-      END IF
+      ENDIF
 
 
@@ -660,5 +660,5 @@
 #endif
       !
-      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
+      !!#LB: NO WHY???? IF( nn_ice == 0 ) THEN
       CALL iom_put( "rho_air"  ,   rhoa )                 ! output air density (kg/m^3) !#LB
       CALL iom_put( "qlw_oce"  ,   zqlw )                 ! output downward longwave heat over the ocean
@@ -674,8 +674,8 @@
       CALL iom_put( 'snowpre', sprecip )                 ! Snow
       CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
-      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
+      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
          CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
-      END IF
+      ENDIF
       !!#LB. ENDIF
       !
@@ -819,8 +819,8 @@
          zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
       CASE( np_humi_dpt )
-         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
+         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
          zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
       CASE( np_humi_rlh )
-         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
+         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
          zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file
       END SELECT
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 11963)
@@ -68,5 +68,5 @@
       INTEGER :: ierr
       !!---------------------------------------------------------------------
-      IF ( l_use_wl ) THEN
+      IF( l_use_wl ) THEN
          ierr = 0
          ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
@@ -76,11 +76,11 @@
          dT_wl(:,:)  = 0._wp
          Hz_wl(:,:)  = Hwl_max
-      END IF
-      IF ( l_use_cs ) THEN
+      ENDIF
+      IF( l_use_cs ) THEN
          ierr = 0
          ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
          IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P0_INIT => allocation of dT_cs failed!' )
          dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
-      END IF
+      ENDIF
    END SUBROUTINE sbcblk_algo_coare3p0_init
 
@@ -190,5 +190,5 @@
       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0'
       !!----------------------------------------------------------------------------------
-      IF ( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl)
+      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl)
 
       l_zt_equal_zu = .FALSE.
@@ -197,16 +197,16 @@
 
       !! Initializations for cool skin and warm layer:
-      IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
 
-      IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
 
-      IF ( l_use_cs .OR. l_use_wl ) THEN
+      IF( l_use_cs .OR. l_use_wl ) THEN
          ALLOCATE ( zsst(jpi,jpj) )
          zsst = T_s ! backing up the bulk SST
          IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
+      ENDIF
 
 
@@ -265,5 +265,5 @@
          dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
          dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-      END IF
+      ENDIF
 
       !! ITERATION BLOCK
@@ -288,8 +288,8 @@
             zeta_t = zt*ztmp0
             zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
-         END IF
+         ENDIF
 
          !! Adjustment the wind at 10m (not needed in the current algo form):
-         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
+         !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
 
          !! Roughness lengthes z0, z0t (z0q = z0t) :
@@ -299,5 +299,5 @@
 
          ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
-         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
+         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LB: some use 1.15 not 1.1 !!!
          z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
 
@@ -315,5 +315,5 @@
             t_zu = t_zt - t_star/vkarmn*ztmp1
             q_zu = q_zt - q_star/vkarmn*ztmp1
-         END IF
+         ENDIF
 
 
@@ -330,5 +330,5 @@
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
 
-         END IF
+         ENDIF
 
          IF( l_use_wl ) THEN
@@ -343,10 +343,10 @@
             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
-         END IF
+         ENDIF
 
          IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-         END IF
+         ENDIF
 
       END DO !DO j_itt = 1, nb_itt
@@ -365,9 +365,9 @@
       IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
 
-      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
-      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
-      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
-
-      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
+      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
+      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
+
+      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
 
    END SUBROUTINE turb_coare3p0
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 11963)
@@ -68,5 +68,5 @@
       INTEGER :: ierr
       !!---------------------------------------------------------------------
-      IF ( l_use_wl ) THEN
+      IF( l_use_wl ) THEN
          ierr = 0
          ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
@@ -76,11 +76,11 @@
          dT_wl(:,:)  = 0._wp
          Hz_wl(:,:)  = Hwl_max
-      END IF
-      IF ( l_use_cs ) THEN
+      ENDIF
+      IF( l_use_cs ) THEN
          ierr = 0
          ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
          IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P6_INIT => allocation of dT_cs failed!' )
          dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
-      END IF
+      ENDIF
    END SUBROUTINE sbcblk_algo_coare3p6_init
 
@@ -190,5 +190,5 @@
       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6'
       !!----------------------------------------------------------------------------------
-      IF ( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl)
+      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl)
 
       l_zt_equal_zu = .FALSE.
@@ -197,16 +197,16 @@
 
       !! Initializations for cool skin and warm layer:
-      IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
 
-      IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
 
-      IF ( l_use_cs .OR. l_use_wl ) THEN
+      IF( l_use_cs .OR. l_use_wl ) THEN
          ALLOCATE ( zsst(jpi,jpj) )
          zsst = T_s ! backing up the bulk SST
          IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
+      ENDIF
 
 
@@ -265,5 +265,5 @@
          dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
          dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-      END IF
+      ENDIF
 
       !! ITERATION BLOCK
@@ -288,8 +288,8 @@
             zeta_t = zt*ztmp0
             zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
-         END IF
+         ENDIF
 
          !! Adjustment the wind at 10m (not needed in the current algo form):
-         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
+         !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
 
          !! Roughness lengthes z0, z0t (z0q = z0t) :
@@ -315,5 +315,5 @@
             t_zu = t_zt - t_star/vkarmn*ztmp1
             q_zu = q_zt - q_star/vkarmn*ztmp1
-         END IF
+         ENDIF
 
 
@@ -330,5 +330,5 @@
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
 
-         END IF
+         ENDIF
 
          IF( l_use_wl ) THEN
@@ -343,10 +343,10 @@
             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
-         END IF
+         ENDIF
 
          IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-         END IF
+         ENDIF
 
       END DO !DO j_itt = 1, nb_itt
@@ -365,9 +365,9 @@
       IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
 
-      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
-      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
-      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
-
-      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
+      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
+      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
+
+      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
 
    END SUBROUTINE turb_coare3p6
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 11963)
@@ -74,5 +74,5 @@
       INTEGER :: ierr
       !!---------------------------------------------------------------------
-      IF ( l_use_wl ) THEN
+      IF( l_use_wl ) THEN
          ierr = 0
          ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
@@ -80,11 +80,11 @@
          dT_wl(:,:)  = 0._wp
          Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars)
-      END IF
-      IF ( l_use_cs ) THEN
+      ENDIF
+      IF( l_use_cs ) THEN
          ierr = 0
          ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
          IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' )
          dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
-      END IF
+      ENDIF
    END SUBROUTINE sbcblk_algo_ecmwf_init
 
@@ -195,5 +195,5 @@
       !!----------------------------------------------------------------------------------
 
-      IF ( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
+      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
 
       l_zt_equal_zu = .FALSE.
@@ -201,16 +201,16 @@
 
       !! Initializations for cool skin and warm layer:
-      IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
 
-      IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
+      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
          &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
 
-      IF ( l_use_cs .OR. l_use_wl ) THEN
+      IF( l_use_cs .OR. l_use_wl ) THEN
          ALLOCATE ( zsst(jpi,jpj) )
          zsst = T_s ! backing up the bulk SST
          IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
-      END IF
+      ENDIF
 
 
@@ -270,5 +270,5 @@
          dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
          dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-      END IF
+      ENDIF
 
 
@@ -293,5 +293,5 @@
          Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
          !! Note: it is slightly different that the L we would get with the usual
-         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
+         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...)
 
          !! Update func_m with new Linv:
@@ -335,5 +335,5 @@
             ztmp1  = LOG(zt/zu) + ztmp2
             q_zu   = q_zt - q_star/vkarmn*ztmp1
-         END IF
+         ENDIF
 
          !! Updating because of updated z0 and z0t and new Linv...
@@ -355,5 +355,5 @@
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
 
-         END IF
+         ENDIF
 
          IF( l_use_wl ) THEN
@@ -366,10 +366,10 @@
             IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
-         END IF
+         ENDIF
 
          IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
-         END IF
+         ENDIF
 
       END DO !DO j_itt = 1, nb_itt
@@ -384,9 +384,9 @@
       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
 
-      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
-      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
-      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
-
-      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
+      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
+      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
+
+      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
 
    END SUBROUTINE turb_ecmwf
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 11963)
@@ -171,5 +171,5 @@
             q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
             q_zu = max(0._wp, q_zu)
-         END IF
+         ENDIF
 
          ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90	(revision 11963)
@@ -541,5 +541,5 @@
             pQns(ji,jj) = zQlat + zQsen + zQlw
 
-            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
+            IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
          END DO
       END DO
@@ -594,6 +594,6 @@
             pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
 
-            IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
-            IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
+            IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
+            IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
 
          END DO
@@ -647,6 +647,6 @@
       pQlat = L_vap(pTs) * zevap
 
-      IF ( PRESENT(pEvap) ) pEvap = - zevap
-      IF ( PRESENT(prhoa) ) prhoa = zrho
+      IF( PRESENT(pEvap) ) pEvap = - zevap
+      IF( PRESENT(prhoa) ) prhoa = zrho
 
    END SUBROUTINE BULK_FORMULA_SCLR
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90	(revision 11963)
@@ -2,20 +2,16 @@
    !!======================================================================
    !!                   ***  MODULE  sbcblk_skin_coare  ***
-   !! Computes:
-   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
-   !!     scheme used at ECMWF
-   !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019)
    !!
-   !!   From routine turb_ecmwf maintained and developed in AeroBulk
-   !!            (https://github.com/brodeau/aerobulk)
+   !!   Module that gathers the cool-skin and warm-layer parameterization used
+   !!   in the COARE family of bulk parameterizations.
    !!
-   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!   Based on the last update for version COARE 3.6 (Fairall et al., 2019)
+   !!
+   !!   Module 'sbcblk_skin_coare' also maintained and developed in AeroBulk (as
+   !!            'mod_skin_coare')
+   !!    (https://github.com/brodeau/aerobulk)   !!
+   !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
-   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
-   !!               based on the cool-skin/warm-layer scheme used at ECMWF
+   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -24,8 +20,8 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
 
-   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
-
-   USE sbcdcy          !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)
-   
+   USE sbcblk_phy      ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc)
+
+   USE sbcdcy          !#LB: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)
+
    USE lib_mpp         ! distribued memory computing library
    USE in_out_manager  ! I/O manager
@@ -38,13 +34,17 @@
 
    !! Cool-skin related parameters:
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
-      &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect
+   !                                                      ! => temperature difference between air-sea interface (z=0)
+   !                                                      !    and right below viscous layer (z=delta)
 
    !! Warm-layer related parameters:
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
-      &                        dT_wl,     &  !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
-      &                        Hz_wl,     &  !: depth of warm-layer [m]
-      &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
-      &                        Tau_ac        !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect
+   !                                                      ! => difference between "almost surface (right below
+   !                                                      !    viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m]
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Qnt_ac !: time integral / accumulated heat stored by the warm layer
+   !                                                      !         Qxdt => [J/m^2] (reset to zero every midnight)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Tau_ac !: time integral / accumulated momentum
+   !                                                      !         Tauxdt => [N.s/m^2] (reset to zero every midnight)
 
    REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp    !: maximum depth of warm layer (adjustable)
@@ -85,5 +85,5 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jc
-      REAL(wp) :: zQabs, zdelta, zfr
+      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus
       !!---------------------------------------------------------------------
       DO jj = 1, jpj
@@ -91,15 +91,23 @@
 
             zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
-            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
-
-            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
-
-            DO jc = 1, 4 ! because implicit in terms of zdelta...
-               zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
+            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt...
+
+            zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K]
+            zqlat = pQlat(ji,jj)
+            zus   = pustar(ji,jj)
+
+
+            zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus )
+
+            DO jc = 1, 4 ! because implicit in terms of zdlt...
+               zfr = MAX( 0.137_wp + 11._wp*zdlt &
+                  &       - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) &
+                  &      , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b)
+               !                  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
                zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
-               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
+               zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus )
             END DO
 
-            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
+            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
 
          END DO
@@ -132,5 +140,5 @@
       REAL(wp) :: zdTwl, zHwl, zQabs, zfr
       REAL(wp) :: zqac, ztac
-      REAL(wp) :: zalpha, zcd1, zcd2, flg
+      REAL(wp) :: zalfa, zcd1, zcd2, flg
       !!---------------------------------------------------------------------
 
@@ -142,5 +150,5 @@
       !! INITIALIZATION:
       zQabs  = 0._wp       ! total heat flux absorped in warm layer
-      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!!
+      zfr    = zfr0        ! initial value of solar flux absorption !#LB: save it and use previous value !!!
 
       IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize!
@@ -148,9 +156,4 @@
       ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ...
 
-      IF (lwp) THEN
-         WRITE(numout,*) ''
-         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime
-      END IF
-      
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -166,44 +169,42 @@
 
             !*****  variables for warm layer  ***
-            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
-
-            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1
-            zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
-
-            
+            zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!)
+
+            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalfa*grav*rho0_w))        !mess-o-constants 1
+            zcd2 = SQRT(2._wp*zalfa*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
+
+
             znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp )   ! 0<rnoon<1. => rnoon*24 = UTC time of local noon
             zmidn = MOD(              znoon-0.5_wp                 , 1._wp )
-            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 
-
-            IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN
+            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight
+
+            IF( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN
                ! Dawn reset to 0!
                l_exit       = .TRUE.
                l_destroy_wl = .TRUE.
-            END IF
-
-            IF ( .NOT. l_exit ) THEN
+            ENDIF
+
+            IF( .NOT. l_exit ) THEN
                !! Initial test on initial guess of absorbed heat flux in warm-layer:
-               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
-                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
-               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!!
-
-               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
+               zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer
+               !                                                        ! => #LB: depends of zfr, which is wild guess... Wrong!!!
+               IF( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
                   ! We have not started to build a WL yet (dT==0) and there's no way it can occur now
                   ! since zQabs <= 0._wp
                   ! => no need to go further
                   l_exit = .TRUE.
-               END IF
-
-            END IF
+               ENDIF
+
+            ENDIF
 
             ! Okay test on updated absorbed flux:
-            !LOLO: remove??? has a strong influence !!!
-            IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
+            !#LB: remove??? has a strong influence !!!
+            IF( (.NOT. l_exit).AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
                l_exit       = .TRUE.
                l_destroy_wl = .TRUE.
-            END IF
-
-
-            IF ( .NOT. l_exit) THEN
+            ENDIF
+
+
+            IF( .NOT. l_exit) THEN
 
                ! Two possibilities at this point:
@@ -217,13 +218,11 @@
                !! some part is useless if Qsw=0 !!!
                DO jl = 1, 5
-                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
-                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
-                  zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)
+                  zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj)
                   zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed
-                  IF ( zqac <= 0._wp ) EXIT
+                  IF( zqac <= 0._wp ) EXIT
                   zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
                END DO
 
-               IF ( zqac <= 0._wp ) THEN
+               IF( zqac <= 0._wp ) THEN
                   l_destroy_wl = .TRUE.
                   l_exit       = .TRUE.
@@ -234,9 +233,9 @@
                   flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl)
                   zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl )
-               END IF
-
-            END IF !IF ( .NOT. l_exit)
-
-            IF ( l_destroy_wl ) THEN
+               ENDIF
+
+            ENDIF !IF( .NOT. l_exit)
+
+            IF( l_destroy_wl ) THEN
                zdTwl = 0._wp
                zfr   = 0.75_wp
@@ -244,7 +243,7 @@
                zqac  = 0._wp
                ztac  = 0._wp
-            END IF
-
-            IF ( iwait == 0 ) THEN
+            ENDIF
+
+            IF( iwait == 0 ) THEN
                !! Iteration loop within bulk algo is over, time to update what needs to be updated:
                dT_wl(ji,jj)  = zdTwl
@@ -252,5 +251,5 @@
                Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
                Tau_ac(ji,jj) = ztac
-            END IF
+            ENDIF
 
          END DO
@@ -276,5 +275,6 @@
       REAL(wp)                :: delta_skin_layer
       REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
-      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
+      REAL(wp), INTENT(in)    :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2]
+      !                              !  => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
       REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
       REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
@@ -282,6 +282,6 @@
       REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
       !!---------------------------------------------------------------------
-      
-      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
+
+      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! #LB: Double check sign + division by palpha !!! units are okay!
 
       ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
@@ -299,4 +299,17 @@
    END FUNCTION delta_skin_layer
 
+
+   FUNCTION frac_solar_abs( pHwl )
+      !!---------------------------------------------------------------------
+      !! Fraction of solar heat flux absorbed inside warm layer
+      !!---------------------------------------------------------------------
+      REAL(wp)             :: frac_solar_abs
+      REAL(wp), INTENT(in) :: pHwl   ! thickness of warm-layer [m]
+      !!---------------------------------------------------------------------
+      frac_solar_abs = 1._wp - ( 0.28*0.014  *(1._wp - EXP(-pHwl/0.014)) &
+         &                       + 0.27*0.357*(1._wp - EXP(-pHwl/0.357)) &
+         &                       + 0.45*12.82*(1-EXP(-pHwl/12.82)) ) / pHwl
+   END FUNCTION frac_solar_abs
+
    !!======================================================================
 END MODULE sbcblk_skin_coare
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90	(revision 11963)
@@ -2,21 +2,31 @@
    !!======================================================================
    !!                   ***  MODULE  sbcblk_skin_ecmwf  ***
-   !! Computes:
-   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
-   !!     scheme used at ECMWF
-   !!    Using formulation/param. of IFS of ECMWF (cycle 40r1)
-   !!   based on IFS doc (avaible online on the ECMWF's website)
-   !!
-   !!   From routine turb_ecmwf maintained and developed in AeroBulk
-   !!            (https://github.com/brodeau/aerobulk)
-   !!
-   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
+   !!
+   !!   Module that gathers the cool-skin and warm-layer parameterization used
+   !!   by the IFS at ECMWF (recoded from scratch =>
+   !!   https://github.com/brodeau/aerobulk)
+   !!
+   !!  Mainly based on Zeng & Beljaars, 2005 with the more recent add-up from
+   !!  Takaya et al., 2010 when it comes to the warm-layer parameterization
+   !!  (contribution of extra mixing due to Langmuir circulation)
+   !!
+   !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin
+   !!    temperature for modeling and data assimilation. Geophysical Research
+   !!    Letters, 32 (14) , pp. 1-4.
+   !!
+   !!  - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen,
+   !!    2010: Refinements to a prognostic scheme of skin sea surface
+   !!    temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985
+   !!
+   !!   Most of the formula are taken from the documentation of IFS of ECMWF
+   !!            (cycle 40r1) (avaible online on the ECMWF's website)
+   !!
+   !!   Routine 'sbcblk_skin_ecmwf' also maintained and developed in AeroBulk (as
+   !!            'mod_skin_ecmwf')
+   !!    (https://github.com/brodeau/aerobulk)
+   !!
+   !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
    !!----------------------------------------------------------------------
-   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
-   !!               based on the cool-skin/warm-layer scheme used at ECMWF
+   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -25,5 +35,5 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
 
-   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
+   USE sbcblk_phy      ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc)
 
    USE lib_mpp         ! distribued memory computing library
@@ -31,5 +41,4 @@
    USE lib_fortran     ! to use key_nosignedzero
 
-
    IMPLICIT NONE
    PRIVATE
@@ -38,11 +47,13 @@
 
    !! Cool-skin related parameters:
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
-      &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect
+   !                                                      ! => temperature difference between air-sea interface (z=0)
+   !                                                      !    and right below viscous layer (z=delta)
 
    !! Warm-layer related parameters:
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
-      &                        dT_wl,  &   !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
-      &                        Hz_wl       !: depth of warm-layer [m]
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect
+   !                                                      ! => difference between "almost surface (right below
+   !                                                      !    viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m]
    !
    REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
@@ -62,8 +73,7 @@
       !! Cool-skin parameterization, based on Fairall et al., 1996:
       !!
-      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
-      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
-      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
-      !! doi:10.1029/95JC03190.
+      !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface
+      !!    skin temperature for modeling and data assimilation. Geophysical
+      !!    Research Letters, 32 (14) , pp. 1-4.
       !!
       !!------------------------------------------------------------------
@@ -81,5 +91,5 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jc
-      REAL(wp) :: zQabs, zdelta, zfr
+      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus
       !!---------------------------------------------------------------------
       DO jj = 1, jpj
@@ -87,16 +97,21 @@
 
             zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
-            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
-
-            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
-
-            DO jc = 1, 4 ! because implicit in terms of zdelta...
-               zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
-               !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
+            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt...
+
+            zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K]
+            zus   = pustar(ji,jj)
+
+            zdlt = delta_skin_layer( zalfa, zQabs, zus )
+
+            DO jc = 1, 4 ! because implicit in terms of zdlt...
+               zfr = MAX( 0.065_wp + 11._wp*zdlt  &
+                  &       - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) &
+                  &      , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
+               !                     !  =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
                zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
-               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
+               zdlt = delta_skin_layer( zalfa, zQabs, zus )
             END DO
 
-            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
+            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
 
          END DO
@@ -106,14 +121,22 @@
 
 
-
    SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST,  pustk )
       !!---------------------------------------------------------------------
       !!
-      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
-      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
+      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) with the
+      !!  more recent add-up from Takaya et al., 2010 when it comes to the
+      !!  warm-layer parameterization (contribution of extra mixing due to
+      !!  Langmuir circulation)
+      !!
+      !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin
+      !!    temperature for modeling and data assimilation. Geophysical Research
+      !!    Letters, 32 (14) , pp. 1-4.
+      !!
+      !!  - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen,
+      !!    2010: Refinements to a prognostic scheme of skin sea surface
+      !!    temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985
       !!
       !!  STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER!
       !!
-      !!    As included in IFS Cy45r1   /  E.C.M.W.F.
       !!     ------------------------------------------------------------------
       !!
@@ -133,15 +156,14 @@
       INTEGER :: ji, jj, jc
       !
-      REAL(wp) :: &
-         & zHwl,    &  !: thickness of the warm-layer [m]
-         & ztcorr,  &  !: correction of dT w.r.t measurement depth of bulk SST (first T-point)
-         & zalpha, & !: thermal expansion coefficient of sea-water [1/K]
-         & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
-         & zfr, zeta, &
-         & zusw, zusw2, &
-         & zLa, zfLa, &
-         & flg, zwf, zQabs, &
-         & ZA, ZB, zL1, zL2, &
-         &  zcst0, zcst1, zcst2, zcst3
+      REAL(wp) :: zHwl      !: thickness of the warm-layer [m]
+      REAL(wp) :: ztcorr    !: correction of dT w.r.t measurement depth of bulk SST (first T-point)
+      REAL(wp) :: zalfa     !: thermal expansion coefficient of sea-water [1/K]
+      REAL(wp) :: zdTwl_b, zdTwl_n  !: temp. diff. between "almost surface (right below viscous layer) and bottom of WL
+      REAL(wp) :: zfr, zeta 
+      REAL(wp) :: zusw, zusw2 
+      REAL(wp) :: zLa, zfLa 
+      REAL(wp) :: flg, zwf, zQabs 
+      REAL(wp) :: ZA, ZB, zL1, zL2
+      REAL(wp) :: zcst0, zcst1, zcst2, zcst3
       !
       LOGICAL :: l_pustk_known
@@ -149,5 +171,5 @@
 
       l_pustk_known = .FALSE.
-      IF ( PRESENT(pustk) ) l_pustk_known = .TRUE.
+      IF( PRESENT(pustk) ) l_pustk_known = .TRUE.
 
       DO jj = 1, jpj
@@ -158,5 +180,5 @@
 
             !! Previous value of dT / warm-layer, adapted to depth:
-            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$
+            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$
             ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl
             zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp )
@@ -166,6 +188,5 @@
             !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl !
 
-            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
-
+            zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!)
 
             ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
@@ -178,22 +199,22 @@
 
             ! Langmuir:
-            IF ( l_pustk_known ) THEN
+            IF( l_pustk_known ) THEN
                zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6))
             ELSE
                zla = 0.3_wp
-            END IF
+            ENDIF
             zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6)
 
             zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0.
 
-            zcst1 = vkarmn*grav*zalpha
+            zcst1 = vkarmn*grav*zalfa
 
             ! 1/L when zQabs > 0 :
             zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw)
-               
+
             zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ???
 
             zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl
-            
+
             ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w )
 
@@ -203,5 +224,5 @@
             !! the sake of optimizations... So all these operations are not done
             !! over and over within the iteration loop...
-            
+
             !! T R U L L Y   I M P L I C I T => needs itteration
             !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0..
@@ -209,11 +230,10 @@
             zdTwl_n = zdTwl_b
             DO jc = 1, 10
-               
+
                zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence
-               
+
                ! 1/L when zdTwl > 0 .AND. zQabs < 0 :
-               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw
-               !zL1 = vkarmn*SQRT( zdTwl_n       *grav*zalpha        / ( 5._wp*zHwl ) ) / zusw   ! => vkarmn outside, not inside zcst1 (just for this particular line) ???
-               
+               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalfa/( 5._wp*zHwl ) ) / zusw
+
                ! Stability parameter (z/L):
                zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2
@@ -224,13 +244,12 @@
 
             END DO
-            
+
             !! Update:
             dT_wl(ji,jj) = zdTwl_n * ztcorr
-            
+
          END DO
       END DO
 
    END SUBROUTINE WL_ECMWF
-
 
 
@@ -249,5 +268,6 @@
       REAL(wp)                :: delta_skin_layer
       REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
-      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
+      REAL(wp), INTENT(in)    :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2]
+      !                              !  => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
       REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
       !!---------------------------------------------------------------------
@@ -255,6 +275,6 @@
       !!---------------------------------------------------------------------
       ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
-      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
-      !
+      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1
+      !                                 !    (ex: weak evaporation and strong positive sensible heat flux)
       zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
       zusw2 = zusw*zusw
@@ -281,12 +301,9 @@
       REAL(wp) :: ztf, zzt2
       !!---------------------------------------------------------------------
-      !
       zzt2 = pzeta*pzeta
-      !
       ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1
       !                                   ! zeta < 0 => ztf = 0
       PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0
          &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0
-      !
    END FUNCTION PHI
 
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbccpl.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbccpl.F90	(revision 11963)
@@ -453,5 +453,5 @@
       CASE( 'conservative'  )
          srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
-         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
+         IF( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
       CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
       CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
@@ -558,5 +558,5 @@
       srcv(jpr_botm )%clname = 'OBotMlt'
       IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
-         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
+         IF( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
             srcv(jpr_topm:jpr_botm)%nct = nn_cats_cpl
          ELSE
@@ -569,7 +569,7 @@
       !                                                      ! ------------------------- !
       srcv(jpr_ts_ice)%clname = 'OTsfIce'    ! needed by Met Office
-      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE.
-      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl
-      IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
+      IF( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE.
+      IF( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl
+      IF( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
 
       !                                                      ! ------------------------- !
@@ -693,5 +693,5 @@
          ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
          DO jn = 1, jprcv
-            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
+            IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
          END DO
          !
@@ -720,18 +720,18 @@
       ! =================================================== !
       DO jn = 1, jprcv
-         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
+         IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
       END DO
       ! Allocate taum part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
+      IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
       ! Allocate w10m part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
+      IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
       ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
-      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
-      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
+      IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
+      IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
       ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
       IF( k_ice /= 0 ) THEN
-         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
-         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
-      END IF
+         IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
+         IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
+      ENDIF
 
       ! ================================ !
@@ -757,5 +757,5 @@
       CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice' )
          ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
-         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = nn_cats_cpl
+         IF( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = nn_cats_cpl
       CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
@@ -777,5 +777,5 @@
       !     1. sending mixed oce-ice albedo or
       !     2. receiving mixed oce-ice solar radiation 
-      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
+      IF( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
          CALL oce_alb( zaos, zacs )
          ! Due to lack of information on nebulosity : mean clear/overcast sky
@@ -796,9 +796,9 @@
          ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used producing atmos-to-ice fluxes (Met Office requirement)
 ! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
-         IF ( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct  = nn_cats_cpl
-         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl
+         IF( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct  = nn_cats_cpl
+         IF( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl
       ENDIF
       
-      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
+      IF(TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
 
       SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
@@ -806,10 +806,10 @@
       CASE( 'ice and snow' ) 
          ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
-         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
+         IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
             ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl
          ENDIF
       CASE ( 'weighted ice and snow' ) 
          ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
-         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl
+         IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
       END SELECT
@@ -828,9 +828,9 @@
          ssnd(jps_a_p)%laction  = .TRUE. 
          ssnd(jps_ht_p)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
+         IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
             ssnd(jps_a_p)%nct  = nn_cats_cpl 
             ssnd(jps_ht_p)%nct = nn_cats_cpl 
          ELSE 
-            IF ( nn_cats_cpl > 1 ) THEN 
+            IF( nn_cats_cpl > 1 ) THEN 
                CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 
             ENDIF 
@@ -839,5 +839,5 @@
          ssnd(jps_a_p)%laction  = .TRUE. 
          ssnd(jps_ht_p)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
+         IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
             ssnd(jps_a_p)%nct  = nn_cats_cpl  
             ssnd(jps_ht_p)%nct = nn_cats_cpl  
@@ -914,8 +914,8 @@
       CASE ( 'ice only' ) 
          ssnd(jps_ttilyr)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN 
+         IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN 
             ssnd(jps_ttilyr)%nct = nn_cats_cpl 
          ELSE 
-            IF ( nn_cats_cpl > 1 ) THEN 
+            IF( nn_cats_cpl > 1 ) THEN 
                CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_ttilyr%cldes if not exchanging category fields' ) 
             ENDIF 
@@ -923,5 +923,5 @@
       CASE ( 'weighted ice' ) 
          ssnd(jps_ttilyr)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl 
+         IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl 
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_ttilyr%cldes;'//sn_snd_ttilyr%cldes ) 
       END SELECT 
@@ -933,8 +933,8 @@
       CASE ( 'ice only' ) 
          ssnd(jps_kice)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 
+         IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 
             ssnd(jps_kice)%nct = nn_cats_cpl 
          ELSE 
-            IF ( nn_cats_cpl > 1 ) THEN 
+            IF( nn_cats_cpl > 1 ) THEN 
                CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 
             ENDIF 
@@ -942,5 +942,5 @@
       CASE ( 'weighted ice' ) 
          ssnd(jps_kice)%laction = .TRUE. 
-         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl 
+         IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl 
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes ) 
       END SELECT 
@@ -1003,5 +1003,5 @@
          ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
          DO jn = 1, jpsnd
-            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
+            IF( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
          END DO
          !
@@ -1030,5 +1030,5 @@
       CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
       
-      IF (ln_usecplmask) THEN 
+      IF(ln_usecplmask) THEN 
          xcplmask(:,:,:) = 0.
          CALL iom_open( 'cplmask', inum )
@@ -1266,5 +1266,5 @@
     
           IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible) 
-      END IF 
+      ENDIF 
       !
       IF( ln_sdw ) THEN  ! Stokes Drift correction activated
@@ -1415,5 +1415,5 @@
          ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
          ELSE                                       ;   zqns(:,:) = 0._wp
-         END IF
+         ENDIF
          ! update qns over the free ocean with:
          IF( nn_components /= jp_iam_opa ) THEN
@@ -1687,5 +1687,5 @@
       ! --- evaporation over ice (kg/m2/s) --- !
       DO jl=1,jpl
-         IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
+         IF(sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl)
          ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF
       ENDDO
@@ -1786,5 +1786,5 @@
       CASE( 'conservative' )     ! the required fields are directly provided
          zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
-         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
+         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
@@ -1795,5 +1795,5 @@
       CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
          zqns_tot(:,:) =  ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
-         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
+         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
@@ -1897,13 +1897,13 @@
 #endif
       ! outputs
-      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
-      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
-      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
-      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
+      IF( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
+      IF( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
+      IF( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
+      IF( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
            &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
-      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
-      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
+      IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
+      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
            &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean)
-      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
+      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
            &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
@@ -1916,5 +1916,5 @@
       CASE( 'conservative' )
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
-         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
             zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
          ELSE
@@ -1928,5 +1928,5 @@
       CASE( 'oce and ice' )
          zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
-         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
             DO jl = 1, jpl
                zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
@@ -1984,5 +1984,5 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
+         IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
             zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
          ELSE
@@ -2062,5 +2062,5 @@
       IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
          
-         IF ( nn_components == jp_iam_opa ) THEN
+         IF( nn_components == jp_iam_opa ) THEN
             ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
          ELSE
@@ -2467,5 +2467,5 @@
       IF( ssnd(jps_ficet)%laction ) THEN 
          CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
-      END IF 
+      ENDIF 
       !                                                      ! ------------------------- ! 
       !                                                      !   Water levels to waves   ! 
@@ -2482,5 +2482,5 @@
          ENDIF  
          CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
-      END IF 
+      ENDIF 
       !
       !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcdcy.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcdcy.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcdcy.F90	(revision 11963)
@@ -94,5 +94,5 @@
             WRITE(numout,*)
          ENDIF
-      END IF
+      ENDIF
 
       ! Setting parameters for each new day:
@@ -121,5 +121,5 @@
                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
                   ztmpm = zupusd - zlousd
-                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
+                  IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
                   !
                ELSE                                         ! day time in two parts
@@ -135,5 +135,5 @@
                   ztmpm = ztmpm1 + ztmpm2
                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
-                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1
+                  IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1
                ENDIF
             ELSE                                   ! 24h light or 24h night
@@ -206,10 +206,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
+               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                   ! When is it night?
                   ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj))
                   ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx )
                   ! is it dawn or dusk?
-                  IF ( ztest > 0._wp ) THEN
+                  IF( ztest > 0._wp ) THEN
                      rdawn_dcy(ji,jj) = ztx
                      rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) )
@@ -232,7 +232,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
+               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
                   rscal(ji,jj) = 0.0_wp
-                  IF ( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part
+                  IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part
                      IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN
                         rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
@@ -247,5 +247,5 @@
                   ENDIF
                ELSE
-                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
+                  IF( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
                      rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
                      rscal(ji,jj) = 1._wp / rscal(ji,jj)
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcice_cice.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcice_cice.F90	(revision 11963)
@@ -132,5 +132,5 @@
          IF      ( ksbc == jp_flx ) THEN
             CALL cice_sbc_force(kt)
-         ELSE IF ( ksbc == jp_purecpl ) THEN
+         ELSE IF( ksbc == jp_purecpl ) THEN
             CALL sbc_cpl_ice_flx( fr_i )
          ENDIF
@@ -140,5 +140,5 @@
          CALL cice_sbc_out ( kt, ksbc )
 
-         IF ( ksbc == jp_purecpl )  CALL cice_sbc_hadgam(kt+1)
+         IF( ksbc == jp_purecpl )  CALL cice_sbc_hadgam(kt+1)
 
       ENDIF                                          ! End sea-ice time step only
@@ -168,5 +168,5 @@
       ! there is no restart file.
       ! Values from a CICE restart file would overwrite this
-      IF ( .NOT. ln_rstart ) THEN    
+      IF( .NOT. ln_rstart ) THEN    
          CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
       ENDIF  
@@ -177,10 +177,10 @@
 
 ! Do some CICE consistency checks
-      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
-         IF ( calc_strair .OR. calc_Tsfc ) THEN
+      IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
+         IF( calc_strair .OR. calc_Tsfc ) THEN
             CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
          ENDIF
-      ELSEIF (ksbc == jp_blk) THEN
-         IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
+      ELSEIF(ksbc == jp_blk) THEN
+         IF( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
             CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
          ENDIF
@@ -202,5 +202,5 @@
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
-      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
+      IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
@@ -297,5 +297,5 @@
 ! forced and coupled case 
 
-      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
+      IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
 
          ztmpn(:,:,:)=0.0
@@ -322,5 +322,5 @@
 
 ! Surface downward latent heat flux (CI_5)
-         IF (ksbc == jp_flx) THEN
+         IF(ksbc == jp_flx) THEN
             DO jl=1,ncat
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
@@ -332,5 +332,5 @@
             DO jj=1,jpj
                DO ji=1,jpi
-                  IF (fr_i(ji,jj).eq.0.0) THEN
+                  IF(fr_i(ji,jj).eq.0.0) THEN
                      DO jl=1,ncat
                         ztmpn(ji,jj,jl)=0.0
@@ -351,5 +351,5 @@
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF(ksbc == jp_flx) THEN
                ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
             ELSE
@@ -360,5 +360,5 @@
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF(ksbc == jp_flx) THEN
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
             ELSE
@@ -368,5 +368,5 @@
          ENDDO
 
-      ELSE IF (ksbc == jp_blk) THEN
+      ELSE IF(ksbc == jp_blk) THEN
 
 ! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
@@ -546,5 +546,5 @@
 ! Freshwater fluxes 
 
-      IF (ksbc == jp_flx) THEN
+      IF(ksbc == jp_flx) THEN
 ! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
 ! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
@@ -552,7 +552,7 @@
 ! Better to use evap and tprecip? (but for now don't read in evap in this case)
          emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      ELSE IF (ksbc == jp_blk) THEN
+      ELSE IF(ksbc == jp_blk) THEN
          emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)        
-      ELSE IF (ksbc == jp_purecpl) THEN
+      ELSE IF(ksbc == jp_purecpl) THEN
 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 
 ! This is currently as required with the coupling fields from the UM atmosphere
@@ -584,10 +584,10 @@
 ! Scale qsr and qns according to ice fraction (bulk formulae only)
 
-      IF (ksbc == jp_blk) THEN
+      IF(ksbc == jp_blk) THEN
          qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
          qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
       ENDIF
 ! Take into account snow melting except for fully coupled when already in qns_tot
-      IF (ksbc == jp_purecpl) THEN
+      IF(ksbc == jp_purecpl) THEN
          qsr(:,:)= qsr_tot(:,:)
          qns(:,:)= qns_tot(:,:)
@@ -624,5 +624,5 @@
 
       CALL cice2nemo(aice,fr_i,'T', 1. )
-      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
+      IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
@@ -879,5 +879,5 @@
 !     B. Gather pn into global array (png)
 
-      IF ( jpnij > 1) THEN
+      IF( jpnij > 1) THEN
          CALL mppsync
          CALL mppgather (pn,0,png) 
@@ -892,5 +892,5 @@
 ! (may be OK but not 100% sure)
 
-      IF (nproc==0) THEN     
+      IF(nproc==0) THEN     
 !        pcg(:,:)=0.0
          DO jn=1,jpnij
@@ -1015,5 +1015,5 @@
 ! the lbclnk call on pn will replace these with sensible values
 
-      IF (nproc==0) THEN
+      IF(nproc==0) THEN
          png(:,:,:)=0.0
          DO jn=1,jpnij
@@ -1028,5 +1028,5 @@
 !     C. Scatter png into NEMO field (pn) for each processor
 
-      IF ( jpnij > 1) THEN
+      IF( jpnij > 1) THEN
          CALL mppsync
          CALL mppscatter (png,0,pn) 
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcisf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcisf.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcisf.F90	(revision 11963)
@@ -303,5 +303,5 @@
       !
       ! Allocate public variable
-      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
+      IF( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
       !
       ! initialisation
@@ -440,5 +440,5 @@
             !! Initialize arrays to 0 (each step)
             zt_sum = 0.e0_wp
-            IF ( ik > 1 ) THEN
+            IF( ik > 1 ) THEN
                ! 1. -----------the average temperature between 200m and 600m ---------------------
                DO jk = misfkt(ji,jj),misfkb(ji,jj)
@@ -459,5 +459,5 @@
             ELSE
                qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
-            END IF
+            ENDIF
          END DO
       END DO
@@ -496,5 +496,5 @@
       ! coeficient for linearisation of potential tfreez
       ! Crude approximation for pressure (but commonly used)
-      IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
+      IF( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
          zlamb1 =-0.0564_wp
          zlamb2 = 0.0773_wp
@@ -558,5 +558,5 @@
                   ! compute s freeze
                   zsfrz=(-zbqe-SQRT(zdis))*zaqer
-                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
+                  IF( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
 
                   ! compute t freeze (eq. 22)
@@ -578,18 +578,18 @@
 
          ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
-         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
+         IF( nn_gammablk <  2 ) THEN ; lit = .FALSE.
          ELSE                           
             ! check total number of iteration
-            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            IF(nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
             ELSE                 ; nit = nit + 1
-            END IF
+            ENDIF
 
             ! compute error between 2 iterations
             ! if needed save gammat and compute zhtflx_b for next iteration
             zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
-            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
+            IF( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
             ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
-            END IF
-         END IF
+            ENDIF
+         ENDIF
       END DO
       !
@@ -718,5 +718,5 @@
                   pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
                   pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
-               END IF
+               ENDIF
             END DO
          END DO
@@ -757,5 +757,5 @@
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbku(ji,jj)
-                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+                  IF( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -789,5 +789,5 @@
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbkv(ji,jj)
-                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -869,5 +869,5 @@
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt, mbkt(ji,jj)
-                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
@@ -879,5 +879,5 @@
             END DO
          END DO
-      END IF 
+      ENDIF 
       !
       !==   ice shelf melting distributed over several levels   ==!
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcmod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcmod.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcmod.F90	(revision 11963)
@@ -243,5 +243,5 @@
          fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp
          fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
-      END IF
+      ENDIF
       IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
          IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
@@ -399,5 +399,5 @@
          emp_b (:,:) = emp (:,:)
          sfx_b (:,:) = sfx (:,:)
-         IF ( ln_rnf ) THEN
+         IF( ln_rnf ) THEN
             rnf_b    (:,:  ) = rnf    (:,:  )
             rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
@@ -437,5 +437,5 @@
       IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
       !
-      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves 
+      IF( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves 
       !
       !                                            !==  Misc. Options  ==!
@@ -471,5 +471,5 @@
 !!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
 !!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. )
-      IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
+      IF( ll_wd ) THEN     ! If near WAD point limit the flux for now
          zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
          zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcrnf.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcrnf.F90	(revision 11963)
@@ -439,5 +439,5 @@
          !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
          !
-         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
+         IF( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
             &                                              'be spread through depth by ln_rnf_depth'               )
          !
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbctide.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbctide.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbctide.F90	(revision 11963)
@@ -72,5 +72,5 @@
          ! Temporarily set nsec_day to beginning of day.
          nsec_day_orig = nsec_day
-         IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN 
+         IF( nsec_day /= NINT(0.5_wp * rdt) ) THEN 
             kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
             nsec_day = NINT(0.5_wp * rdt)
Index: /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/tideini.F90
===================================================================
--- /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/tideini.F90	(revision 11962)
+++ /NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/tideini.F90	(revision 11963)
@@ -68,5 +68,5 @@
       !
       IF( ln_tide ) THEN
-         IF (lwp) THEN
+         IF(lwp) THEN
             WRITE(numout,*)
             WRITE(numout,*) 'tide_init : Initialization of the tidal components'
@@ -127,5 +127,5 @@
       kt_tide = nit000
       !
-      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
+      IF(.NOT.ln_scal_load ) rn_scal_load = 0._wp
       !
    END SUBROUTINE tide_init
