Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11858)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11873)
@@ -45,5 +45,5 @@
               &          , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice  & 
               &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      & 
-			  &          , ptaui_ice,     ptauj_ice                      &
+              &          , ptaui_ice,     ptauj_ice                      &
 #endif			  
 			  &   )
@@ -127,11 +127,11 @@
       DO jj = 1,jpj
          DO ji = 1,jpi
-		    zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
-#if defined key_si3				
-			zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
-		    ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
+            zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
+#if defined key_si3
+            zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
+            ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
 #else
-	        ustar2(ji,jj) = zzoce	
-#endif						
+            ustar2(ji,jj) = zzoce	
+#endif
          END DO
       END DO  
@@ -153,7 +153,7 @@
          DO jk = 3, jpkam1
             DO ji = 1, jpi   ! vector opt.
-               z_elem_a( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal
-               z_elem_c( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal       
-               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal           
+               z_elem_a( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
+               z_elem_c( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal       
+               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal           
             END DO
          END DO  
@@ -168,9 +168,9 @@
             z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
          END DO            
-                   
+
          DO jtra = 1,jptq  ! loop on active tracers
                
             DO jk = 3, jpkam1
-               DO ji = fs_2, fs_jpim1                 
+               DO ji = 1,jpi
                   tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
                END DO
@@ -178,25 +178,25 @@
 
             IF(jtra == jp_ta) THEN
-               DO ji = 1,jpi  ! surface boundary condition for temperature              
+               DO ji = 1,jpi  ! boundary conditions for temperature              
                   zztmp1 = psen(ji, jj) 
-				  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )	
+                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )	
 #if defined key_si3				  
                   zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
                   zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
 #endif				  
-				  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1				  
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2               
-				  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
+                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1				  
+                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rdt_abl * zztmp2               
+                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
                END DO 
             ELSE   
-               DO ji = 1,jpi  ! surface boundary condition for humidity              
+               DO ji = 1,jpi  ! boundary conditions for humidity              
                   zztmp1 = pevp(ji, jj) 
-				  zztmp2 = pevp(ji, jj) * pssq(ji, jj)   
+                  zztmp2 = pevp(ji, jj) * pssq(ji, jj)   
 #if defined key_si3				  
                   zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 
-				  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)    
+                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)    
 #endif	
-                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1
-                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2               
+                  z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1
+                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rdt_abl * zztmp2               
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
                END DO                                     
@@ -442,6 +442,6 @@
                   zcff  = (1._wp-zmsk) + zmsk * zcff2 * rdt   ! zcff = 1 for masked points
                                                               ! rdt = rdt_abl / nn_fsbc                          
-		  zcff  = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff
-				  
+                  zcff  = zcff * rest_eq(ji,jj)
+                  z_cft( ji, jj, jk ) = zcff
                   u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
                      &                               + zcff   * pu_dta( ji, jj, jk       )                      
@@ -468,4 +468,5 @@
                zcff  = (1._wp-zmsk) + zmsk * zcff2 * rdt   ! zcff = 1 for masked points
                                                            ! rdt = rdt_abl / nn_fsbc                          
+               !z_cft( ji, jj, jk ) = zcff
                tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
                   &                                       + zcff   * pt_dta( ji, jj, jk )
@@ -474,5 +475,4 @@
                   &                                       + zcff   * pq_dta( ji, jj, jk )
                
-               ! z_cft( ji, jj, jk ) = zcff
             END DO
          END DO
@@ -483,14 +483,40 @@
       !                            !  6 *** MPI exchanges  
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-	  !
-      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1., kfillmode = jpfillnothing  )
-      CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T',  1., tq_abl(:,:,:,nt_a,jp_qa), 'T',  1., kfillmode = jpfillnothing  )   ! ++++ this should not be needed...
-	  !
-      CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      ))
-      CALL iom_put ( 'v_abl',  v_abl(:,:,2:jpka,nt_a      ))
-      CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
-      CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
-      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           ))
-      CALL iom_put (  'pblh',   pblh(:,:                  ))
+      !
+      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1. )
+      !
+      ! first ABL level
+      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) )
+      IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) )
+      IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl",  tq_abl(:,:,2,nt_a,jp_ta) )
+      IF ( iom_use("qz1_abl") ) CALL iom_put ( "qz1_abl",  tq_abl(:,:,2,nt_a,jp_qa) )
+      IF ( iom_use("uz1_dta") ) CALL iom_put ( "uz1_dta",  pu_dta(:,:,2           ) )
+      IF ( iom_use("vz1_dta") ) CALL iom_put ( "vz1_dta",  pv_dta(:,:,2           ) )
+      IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta",  pt_dta(:,:,2           ) )
+      IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta",  pq_dta(:,:,2           ) )
+      ! all ABL levels
+      IF ( iom_use("u_abl"  ) ) CALL iom_put ( "u_abl"  ,   u_abl(:,:,2:jpka,nt_a      ) )
+      IF ( iom_use("v_abl"  ) ) CALL iom_put ( "v_abl"  ,   v_abl(:,:,2:jpka,nt_a      ) )
+      IF ( iom_use("t_abl"  ) ) CALL iom_put ( "t_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_ta) )
+      IF ( iom_use("q_abl"  ) ) CALL iom_put ( "q_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_qa) )
+      IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a      ) )
+      IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) )
+      IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka           ) )
+      IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka           ) )
+      IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) )
+      ! debug (to be removed)
+      IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta",  pu_dta(:,:,2:jpka) )
+      IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta",  pv_dta(:,:,2:jpka) )
+      IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) )
+      IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) )
+      IF ( iom_use("coeft") ) CALL iom_put ( "coeft",   z_cft(:,:,2:jpka) )
+      IF( ln_geos_winds ) THEN
+         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) )
+         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2           ) )
+      END IF
+      IF( ln_hpgls_frc ) THEN
+         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   )
+         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   )
+      END IF
       !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -772,14 +798,14 @@
             END IF 
          END DO
-         ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
-         !
       !-------------      
       END DO      
       !-------------
-	  IF(ln_smth_pblh) THEN
-	     CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
-	     CALL smooth_pblh( pblh, msk_abl )
-	     CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)	  
-	  ENDIF
+      !
+      ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
+      IF(ln_smth_pblh) THEN
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)
+         CALL smooth_pblh( pblh, msk_abl )
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)	  
+      ENDIF
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  Diagnostic mixing length computation
