Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diawri.F90	(revision 12150)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diawri.F90	(revision 12182)
@@ -28,4 +28,5 @@
    USE isf_oce
    USE isfcpl
+   USE abl            ! abl variables in case ln_abl = .true.
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -68,5 +69,7 @@
    PUBLIC   dia_wri_state
    PUBLIC   dia_wri_alloc           ! Called by nemogcm module
-
+#if ! defined key_iomput   
+   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.)
+#endif
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
    INTEGER ::          nb_T              , ndim_bT   ! grid_T file
@@ -74,6 +77,8 @@
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
    INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
+   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file   
    INTEGER ::   ndex(1)                              ! ???
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
@@ -417,8 +422,15 @@
          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
          !
-      dia_wri_alloc = MAXVAL(ierr)
+	  dia_wri_alloc = MAXVAL(ierr)
       CALL mpp_sum( 'diawri', dia_wri_alloc )
       !
    END FUNCTION dia_wri_alloc
+ 
+   INTEGER FUNCTION dia_wri_alloc_abl()
+      !!----------------------------------------------------------------------
+	  ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl)
+      CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
+      !
+   END FUNCTION dia_wri_alloc_abl
 
    
@@ -444,4 +456,5 @@
       INTEGER  ::   ierr                                     ! error code return from allocation
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      INTEGER  ::   ipka                                     ! ABL
       INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
@@ -449,4 +462,5 @@
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -482,4 +496,5 @@
       ijmi = 1      ;      ijma = jpj
       ipk = jpk
+      IF(ln_abl) ipka = jpkam1
 
       ! define time axis
@@ -584,4 +599,20 @@
             &          "m", ipk, gdepw_1d, nz_W, "down" )
 
+         IF( ln_abl ) THEN 
+         ! Define the ABL grid FILE ( nid_A )
+            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
+            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
+               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+               &          nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
+            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
+               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
+            !                                                            ! Index of ocean points
+			ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
+			zw3d_abl(:,:,:) = 1._wp 
+			CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
+            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
+			DEALLOCATE(zw3d_abl)
+         ENDIF
 
          ! Declare all the output fields as NETCDF variables
@@ -633,5 +664,29 @@
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!
+         !
+         IF( ln_abl ) THEN
+            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )		 			   
+#if defined key_si3
+            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
+#endif
+            CALL histend( nid_A, snc4chunks=snc4set )
+         ENDIF
+         !
          IF( ln_icebergs ) THEN
             CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
@@ -791,5 +846,28 @@
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
-!
+      !
+      IF( ln_abl ) THEN 
+         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
+         IF( ln_mskland )   THEN 
+            DO jk=1,jpka
+               zw3d_abl(:,:,jk) = tmask(:,:,1)
+            END DO       
+         ELSE
+            zw3d_abl(:,:,:) = 1._wp		 
+         ENDIF			
+         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh 
+         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
+         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
+         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
+         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl		 
+         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
+         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
+         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl	
+#if defined key_si3
+         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
+#endif
+         DEALLOCATE(zw3d_abl)
+      ENDIF
+      !
       IF( ln_icebergs ) THEN
          !
@@ -861,4 +939,5 @@
          CALL histclo( nid_V )
          CALL histclo( nid_W )
+         IF(ln_abl) CALL histclo( nid_A )
       ENDIF
       !
@@ -883,5 +962,4 @@
       INTEGER           , INTENT( in ) ::   Kmm              ! time level index
       CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zisfdebug
       !!
       INTEGER :: inum, jk
@@ -955,4 +1033,10 @@
          CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
       ENDIF
+      IF ( ln_abl ) THEN
+         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind
+         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind
+         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature
+         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
+      ENDIF
  
 #if defined key_si3
