modi_glt_thermo.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!GLT_LIC The GELATO model is a seaice model used in stand-alone or embedded mode. 
!GLT_LIC  It has been developed by Meteo-France. The holder of GELATO is Meteo-France.
!GLT_LIC  
!GLT_LIC  This software is governed by the CeCILL-C license under French law and biding
!GLT_LIC  by the rules of distribution of free software. See the CeCILL-C_V1-en.txt
!GLT_LIC  (English) and CeCILL-C_V1-fr.txt (French) for details. The CeCILL is a free
!GLT_LIC  software license, explicitly compatible with the GNU GPL
!GLT_LIC  (see http://www.gnu.org/licenses/license-list.en.html#CeCILL)
!GLT_LIC  
!GLT_LIC  The CeCILL-C licence agreement grants users the right to modify and re-use the
!GLT_LIC  software governed by this free software license. The exercising of this right
!GLT_LIC  is conditional upon the obligation to make available to the community the
!GLT_LIC  modifications made to the source code of the software so as to contribute to
!GLT_LIC  its evolution.
!GLT_LIC  
!GLT_LIC  In consideration of access to the source code and the rights to copy, modify
!GLT_LIC  and redistribute granted by the license, users are provided only with a limited
!GLT_LIC  warranty and the software's author, the holder of the economic rights, and the
!GLT_LIC  successive licensors only have limited liability. In this respect, the risks
!GLT_LIC  associated with loading, using, modifying and/or developing or reproducing the
!GLT_LIC  software by the user are brought to the user's attention, given its Free
!GLT_LIC  Software status, which may make it complicated to use, with the result that its
!GLT_LIC  use is reserved for developers and experienced professionals having in-depth
!GLT_LIC  computer knowledge. Users are therefore encouraged to load and test the
!GLT_LIC  suitability of the software as regards their requirements in conditions enabling
!GLT_LIC  the security of their systems and/or data to be ensured and, more generally, to
!GLT_LIC  use and operate it in the same conditions of security. 
!GLT_LIC  
!GLT_LIC  The GELATO sofware is cureently distibuted with the SURFEX software, available at 
!GLT_LIC  http://www.cnrm.meteo.fr/surfex. The fact that you download the software deemed that
!GLT_LIC  you had knowledge of the CeCILL-C license and that you accept its terms.
!GLT_LIC  Attempts to use this software in a way not complying with CeCILL-C license
!GLT_LIC  may lead to prosecution. 
!GLT_LIC 
! =======================================================================
! ======================= MODULE modi_glt_thermo ========================
! =======================================================================
!
!
! Goal:
! -----
!   Computation of the thermodynamic forcing over open water, sea ice
! and snow-covered sea ice. Version with thickness within each box
! which is constant, compensated with variable fractional ice
! covers.
!   This part of the code runs on a reduced grid, the work grid only includes 
! ocean points where SST is less than a certain threshold or where there
! is sea ice. We consider that:
!       - if there is a lot of sea ice, there is a possibility for warm 
! water advection in this area, hence any place where there is sea ice 
! should be treated to avoid losing sea ice
!       - a reasonable temperature threshold is 1°, to include areas 
! where salinity is zero (freezing point = 0°C)
!  
! Created : 2009/06 - Reduced grid option introduced - full grid still 
!   available (D. Salas y Melia)
! Modified: 2009/11 - Full grid is no longer treated (D. Salas y Melia)
! Modified: 2010/05 - mpi Multi-processing introduced (S.Senesi)
! Modified: 2011/04 - (Re)-ensure portability on both NEC and PC
! Modified: 2011/05 - Solve issue occurring when less ice points than 
!   processors (S.Senesi)
! Modified: 2011/12 - Introduce new diagnostics (A. Voldoire)
! Modified: 2012/07 (D. Salas y Melia)
!           Input now consists of fields defined on subdomains
! Modified: 2012/11 (D. Salas y Melia)
!           Sea ice damping
! Modified: 2015/07 (D. Salas y Melia)
!           Sea ice damping has been externalized from thermo_r routine
!
! -------------------- BEGIN MODULE modi_glt_thermo -------------------------

!THXS_SFX!MODULE modi_glt_thermo
!THXS_SFX!INTERFACE
!THXS_SFX!
!THXS_SFX!SUBROUTINE glt_thermo  &
!THXS_SFX!  ( tpdom,pustar,tpmxl,tpatm,  &
!THXS_SFX!    tpblkw,tpblki,tpbud,tpdia,tptfl,tpsit,tpsil,tpsit_d )
!THXS_SFX!  USE modd_types_glt
!THXS_SFX!  USE modd_glt_param
!THXS_SFX!  TYPE(t_dom), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpdom
!THXS_SFX!  REAL, DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        pustar
!THXS_SFX!  TYPE(t_mxl), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpmxl
!THXS_SFX!  TYPE(t_atm), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpatm
!THXS_SFX!  TYPE(t_blk), DIMENSION(nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpblkw
!THXS_SFX!  TYPE(t_blk), DIMENSION(nt,nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpblki
!THXS_SFX!  TYPE(t_bud), DIMENSION(nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpbud
!THXS_SFX!  TYPE(t_dia), DIMENSION(nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpdia
!THXS_SFX!  TYPE(t_tfl), DIMENSION(nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tptfl
!THXS_SFX!  TYPE(t_sit), DIMENSION(nt,nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpsit
!THXS_SFX!  TYPE(t_vtp), DIMENSION(nl,nt,nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpsil
!THXS_SFX!  TYPE(t_sit), DIMENSION(ntd,nx,ny), OPTIONAL, INTENT(in) ::  &
!THXS_SFX!        tpsit_d
!THXS_SFX!END SUBROUTINE glt_thermo
!THXS_SFX!
!THXS_SFX!END INTERFACE
!THXS_SFX!END MODULE modi_glt_thermo

! ---------------------- END MODULE modi_glt_thermo -------------------------
!
!
! -----------------------------------------------------------------------
! ----------------------- SUBROUTINE glt_thermo -----------------------------
!
SUBROUTINE glt_thermo  &
  ( tpdom,pustar,tpmxl,tpatm,  &
    tpblkw,tpblki,tpbud,tpdia,tptfl,tpsit,tpsil,tpsit_d )
!
!
! 1. DECLARATIONS
! ===============
!
! 1.1. Module declarations
! ------------------------
!
  USE modd_glt_const_thm
  USE modd_types_glt
  USE modd_glt_param
  USE mode_glt_stats
  USE mode_glt_stats_r
  USE modi_glt_thermo_r
  USE modi_glt_constrain_r
!
  IMPLICIT NONE
!
!
! 1.2. Dummy arguments declarations
! ---------------------------------
!
  TYPE(t_dom), DIMENSION(nx,ny), INTENT(in) ::  &
        tpdom
  REAL, DIMENSION(nx,ny), INTENT(in) ::  &
        pustar
  TYPE(t_mxl), DIMENSION(nx,ny), INTENT(in) ::  &
        tpmxl
  TYPE(t_atm), DIMENSION(nx,ny), INTENT(in) ::  &
        tpatm
  TYPE(t_blk), DIMENSION(nx,ny), INTENT(inout) ::  &
        tpblkw
  TYPE(t_blk), DIMENSION(nt,nx,ny), INTENT(in) ::  &
        tpblki
  TYPE(t_bud), DIMENSION(nx,ny), INTENT(inout) ::  &
        tpbud
  TYPE(t_dia), DIMENSION(nx,ny), INTENT(inout) ::  &
        tpdia
  TYPE(t_tfl), DIMENSION(nx,ny), INTENT(inout) ::  &
        tptfl
  TYPE(t_sit), DIMENSION(nt,nx,ny), INTENT(inout) ::  &
        tpsit
  TYPE(t_vtp), DIMENSION(nl,nt,nx,ny), INTENT(inout) ::  &
        tpsil
  TYPE(t_sit), DIMENSION(ntd,nx,ny), OPTIONAL, INTENT(in) ::  &
        tpsit_d
!
!
! 1.3. Local variables declarations
! ---------------------------------
!
  INTEGER ::  &
        jk,jl
  LOGICAL, DIMENSION(nx,ny) ::  &
        gsel
  INTEGER, DIMENSION(nx,ny) ::  &
        isel
  REAL, DIMENSION(nx,ny) ::  &
        zfsit,zsnflx
  TYPE(t_dia), DIMENSION(nx,ny) ::  &
        tzdia0
  TYPE(t_tfl), DIMENSION(nx,ny) ::  &
        tztfl0
  TYPE(t_sit), DIMENSION(nt,nx,ny) ::  &
        tzsit0
  TYPE(t_vtp), DIMENSION(nl,nt,nx,ny) ::  &
        tzsil0
!  INTEGER ::  &
!        ji,jj
!  INTEGER, DIMENSION(2,nx,ny) ::  &
!        ind
!  INTEGER, DIMENSION(:,:), ALLOCATABLE ::  &
!        ind_r
  TYPE(t_dom), DIMENSION(:), ALLOCATABLE ::  &
        tzdom_r
  REAL, DIMENSION(:), ALLOCATABLE ::  &
        zustar_r
  TYPE(t_mxl), DIMENSION(:), ALLOCATABLE ::  &
        tzmxl_r
  TYPE(t_atm), DIMENSION(:), ALLOCATABLE ::  &
        tzatm_r
  TYPE(t_blk), DIMENSION(:), ALLOCATABLE ::  &
        tzblkw_r
  TYPE(t_blk), DIMENSION(:,:), ALLOCATABLE ::  &
        tzblki_r
  TYPE(t_bud), DIMENSION(:), ALLOCATABLE ::  &
        tzbud_r
  TYPE(t_dia), DIMENSION(:), ALLOCATABLE ::  &
        tzdia_r
  TYPE(t_tfl), DIMENSION(:), ALLOCATABLE ::  &
        tztfl_r
  TYPE(t_sit), DIMENSION(:,:), ALLOCATABLE ::  &
        tzsit_d_r
  TYPE(t_sit), DIMENSION(:,:), ALLOCATABLE ::  &
        tzsit_r
  TYPE(t_vtp), DIMENSION(:,:,:), ALLOCATABLE ::  &
        tzsil_r
!
!
! 1.4. Welcome message
! --------------------
!
  IF (lp1) THEN
    WRITE(noutlu,*) ' '
    WRITE(noutlu,*) ' *** LEVEL 3 - SUBROUTINE THERMO'
    WRITE(noutlu,*) ' '
  ENDIF
!
!
!
! 2. THERMODYNAMICS ON REDUCED GRID
! ==================================
!
! 2.1. Define grid point selection criterion
! -------------------------------------------
!
  zfsit(:,:) = glt_iceconcm(tpdom,tpsit) 
  isel(:,:) = 0
!
! In principe, calculations on tpdom%imk==1 would suffice. But it
! would require bounding in the end (thermodynamics would not be 
! purely 1d...)
!
  WHERE( tpdom(:,:)%tmk==1 .AND.  &
  (tpmxl(:,:)%tml <= 1. .OR. zfsit(:,:)>epsil1) )
      isel(:,:) = 1
  ENDWHERE
!
  gsel(:,:) = ( isel(:,:)==1 )
! 
!
! 2.2. Pack all fields and allocate reduced grid arrays
! ------------------------------------------------------
!
! .. Reduced grid dimension
!
  np = sum(isel)
!
  IF (lp3) THEN
    WRITE(noutlu,*)  &
      '**********************************************************'
    WRITE(noutlu,*) 'REDUCED GRID:'
    WRITE(noutlu,1000) gelato_myrank, np, nx*ny, nx, ny
    WRITE(noutlu,*)  &
      '**********************************************************'
  ENDIF
!
! .. Note that on some platforms, every locally-allocated array needs to be
! deallocated explicitly (this is done at the end of this subroutine).
! More over, the deallocations need be done in reverse order.
!
  IF ( np /= 0 ) THEN
!    ALLOCATE( ind_r(2,np) )
    ALLOCATE( tzdom_r(np) )
    ALLOCATE( zustar_r(np) )
    ALLOCATE( tzmxl_r(np) )
    ALLOCATE( tzatm_r(np) )
    ALLOCATE( tzblkw_r(np) )
    ALLOCATE( tzblki_r(nt,np) )
    ALLOCATE( tzbud_r(np) )
    ALLOCATE( tzdia_r(np) )
    ALLOCATE( tztfl_r(np) )
! Note that the error check on ntd dimension will be done in glt_constrain_r
  IF ( present(tpsit_d) ) ALLOCATE( tzsit_d_r(ntd,np) )
    ALLOCATE( tzsit_r(nt,np) )
    ALLOCATE( tzsil_r(nl,nt,np) )
!
!
! 2.3. From global to reduced grid arrays
! ----------------------------------------
!
! Ces tableaux devraient passer en variable globale, pour des recherches
! d'indices...
!ind_r(1,:) = PACK( ind(1,:,:),gsel )
!ind_r(2,:) = PACK( ind(2,:,:),gsel )
!if (np>=6832) write(noutlu,*)'index 6832 =',ind_r(1,6832),ind_r(2,6832)
! .. Domain
!
     tzdom_r%tmk = pack( tpdom%tmk ,gsel )
     tzdom_r%umk = pack( tpdom%umk ,gsel )
     tzdom_r%vmk = pack( tpdom%vmk ,gsel )
     tzdom_r%imk = pack( tpdom%imk ,gsel )
     tzdom_r%indi = pack( tpdom%indi ,gsel )
     tzdom_r%indj = pack( tpdom%indj ,gsel )
     tzdom_r%lon = pack( tpdom%lon ,gsel )
     tzdom_r%lat = pack( tpdom%lat ,gsel )
     tzdom_r%dxc = pack( tpdom%dxc ,gsel )
     tzdom_r%dyc = pack( tpdom%dyc ,gsel )
     tzdom_r%srf = pack( tpdom%srf ,gsel )
!
! .. Friction velocity
!
     zustar_r = pack( pustar, gsel )
!
! .. Mixed layer conditions
!
     tzmxl_r%hco = pack( tpmxl%hco ,gsel )
     tzmxl_r%qoc = pack( tpmxl%qoc ,gsel )
     tzmxl_r%sml = pack( tpmxl%sml ,gsel )
     tzmxl_r%ssh = pack( tpmxl%ssh ,gsel )
     tzmxl_r%tml = pack( tpmxl%tml ,gsel )
     tzmxl_r%mlf = pack( tpmxl%mlf ,gsel )
     tzmxl_r%uml = pack( tpmxl%uml ,gsel )
     tzmxl_r%vml = pack( tpmxl%vml ,gsel )
     tzmxl_r%qml = pack( tpmxl%qml ,gsel )
!
! .. Atmospheric conditions
!
     tzatm_r%lip = pack( tpatm%lip ,gsel )
     tzatm_r%sop = pack( tpatm%sop ,gsel )
     tzatm_r%ztx = pack( tpatm%ztx ,gsel )
     tzatm_r%mty = pack( tpatm%mty ,gsel )
!
! .. Surface fluxes (over ocean)
!     
     tzblkw_r%swa = pack( tpblkw%swa ,gsel )
     tzblkw_r%nsf = pack( tpblkw%nsf ,gsel )
     tzblkw_r%dfl = pack( tpblkw%dfl ,gsel )
     tzblkw_r%eva = pack( tpblkw%eva ,gsel )
!
! .. Surface fluxes (over ice)
!
    DO jk=1,nt
        tzblki_r(jk,:)%swa = pack( tpblki(jk,:,:)%swa ,gsel )
        tzblki_r(jk,:)%nsf = pack( tpblki(jk,:,:)%nsf ,gsel )
        tzblki_r(jk,:)%dfl = pack( tpblki(jk,:,:)%dfl ,gsel )
        tzblki_r(jk,:)%eva = pack( tpblki(jk,:,:)%eva ,gsel )
    END DO
!
! .. Energy budget
!
     tzbud_r%eni = pack( tpbud%eni ,gsel )
     tzbud_r%enn = pack( tpbud%enn ,gsel )
     tzbud_r%bii = pack( tpbud%bii ,gsel )
     tzbud_r%nii = pack( tpbud%nii ,gsel )
     tzbud_r%nli = pack( tpbud%nli ,gsel )
     tzbud_r%hii = pack( tpbud%hii ,gsel )
     tzbud_r%hli = pack( tpbud%hli ,gsel )
     tzbud_r%hio = pack( tpbud%hio ,gsel )
     tzbud_r%hlo = pack( tpbud%hlo ,gsel )
     tzbud_r%wii = pack( tpbud%wii ,gsel )
     tzbud_r%wli = pack( tpbud%wli ,gsel )
     tzbud_r%fwi = pack( tpbud%fwi ,gsel )
     tzbud_r%fwn = pack( tpbud%fwn ,gsel )
     tzbud_r%isi = pack( tpbud%isi ,gsel )
     tzbud_r%isn = pack( tpbud%isn ,gsel )
!
! .. Diagnostics
!
     tzdia_r%uvl = pack( tpdia%uvl ,gsel )
     tzdia_r%vvl = pack( tpdia%vvl ,gsel )
     tzdia_r%aiw = pack( tpdia%aiw ,gsel )
     tzdia_r%asi = pack( tpdia%asi ,gsel )
     tzdia_r%amp = pack( tpdia%amp ,gsel )
     tzdia_r%asn = pack( tpdia%asn ,gsel )
     tzdia_r%cgl = pack( tpdia%cgl ,gsel )
     tzdia_r%dsa = pack( tpdia%dsa ,gsel )
     tzdia_r%dsn = pack( tpdia%dsn ,gsel )
     tzdia_r%dsi = pack( tpdia%dsi ,gsel )
     tzdia_r%dwi = pack( tpdia%dwi ,gsel )
     tzdia_r%lip = pack( tpdia%lip ,gsel )
     tzdia_r%lsi = pack( tpdia%lsi ,gsel )
     tzdia_r%mrb = pack( tpdia%mrb ,gsel )
     tzdia_r%mrt = pack( tpdia%mrt ,gsel )
     tzdia_r%mrl = pack( tpdia%mrl ,gsel )
     tzdia_r%sie = pack( tpdia%sie ,gsel )
     tzdia_r%sne = pack( tpdia%sne ,gsel )
     tzdia_r%sni = pack( tpdia%sni ,gsel )
     tzdia_r%snm = pack( tpdia%snm ,gsel )
     tzdia_r%snml = pack( tpdia%snml ,gsel )
     tzdia_r%sop = pack( tpdia%sop ,gsel )
     tzdia_r%the = pack( tpdia%the ,gsel )
     tzdia_r%tin = pack( tpdia%tin ,gsel )
     tzdia_r%tiw = pack( tpdia%tiw ,gsel )
     tzdia_r%qoi = pack( tpdia%qoi ,gsel )
     tzdia_r%xtr = pack( tpdia%xtr ,gsel )
     tzdia_r%ytr = pack( tpdia%ytr ,gsel )
     tzdia_r%sp1 = pack( tpdia%sp1 ,gsel )
     tzdia_r%sp2 = pack( tpdia%sp2 ,gsel )
     tzdia_r%sui = pack( tpdia%sui ,gsel )
     tzdia_r%sut = pack( tpdia%sut ,gsel )
     tzdia_r%sus = pack( tpdia%sus ,gsel )
     tzdia_r%suw = pack( tpdia%suw ,gsel )
     tzdia_r%s_pr = pack( tpdia%s_pr ,gsel )
     tzdia_r%s_prsn = pack( tpdia%s_prsn ,gsel )
     tzdia_r%o_pr = pack( tpdia%o_pr ,gsel )
     tzdia_r%o_prsn = pack( tpdia%o_prsn ,gsel )
     tzdia_r%subcio = pack( tpdia%subcio ,gsel )
     tzdia_r%snicio = pack( tpdia%snicio ,gsel )
     tzdia_r%hsicio = pack( tpdia%hsicio ,gsel )
     tzdia_r%lmlcio = pack( tpdia%lmlcio ,gsel )
     tzdia_r%salcio = pack( tpdia%salcio ,gsel )
     tzdia_r%l_pr = pack( tpdia%l_pr ,gsel )
     tzdia_r%l_prsn = pack( tpdia%l_prsn ,gsel )
     tzdia_r%sul = pack( tpdia%sul ,gsel )
! These ones should not be packed: unused in thermodynamics, and contain
! data outside the ice domain (the data will be lost if packing):
!   - tzdia_r%ifw
!   - tzdia_r%sic
!   - tzdia_r%sit
!   - tzdia_r%atx,%aty
!   - tzdia_r%otx,%oty
!
! .. Fluxes transmitted to the mixed layer
!
     tztfl_r%lio = pack( tptfl%lio ,gsel )
     tztfl_r%llo = pack( tptfl%llo ,gsel )
     tztfl_r%tio = pack( tptfl%tio ,gsel )
     tztfl_r%tlo = pack( tptfl%tlo ,gsel )
     tztfl_r%sio = pack( tptfl%sio ,gsel )
     tztfl_r%cio = pack( tptfl%cio ,gsel )
     tztfl_r%wio = pack( tptfl%wio ,gsel )
     tztfl_r%wlo = pack( tptfl%wlo ,gsel )
     tztfl_r%xio = pack( tptfl%xio ,gsel )
     tztfl_r%yio = pack( tptfl%yio ,gsel )
!
! .. Damping / restoring data
!
  IF ( present(tpsit_d) ) THEN
    DO jk=1,ntd
      tzsit_d_r(jk,:)%esi = pack( tpsit_d(jk,:,:)%esi, gsel )
      tzsit_d_r(jk,:)%asn = pack( tpsit_d(jk,:,:)%asn, gsel )
      tzsit_d_r(jk,:)%fsi = pack( tpsit_d(jk,:,:)%fsi, gsel )
      tzsit_d_r(jk,:)%hsi = pack( tpsit_d(jk,:,:)%hsi, gsel )
      tzsit_d_r(jk,:)%hsn = pack( tpsit_d(jk,:,:)%hsn, gsel )
      tzsit_d_r(jk,:)%rsn = pack( tpsit_d(jk,:,:)%rsn, gsel )
      tzsit_d_r(jk,:)%tsf = pack( tpsit_d(jk,:,:)%tsf, gsel )
      tzsit_d_r(jk,:)%ssi = pack( tpsit_d(jk,:,:)%ssi, gsel )
      tzsit_d_r(jk,:)%age = pack( tpsit_d(jk,:,:)%age, gsel )
      tzsit_d_r(jk,:)%vmp = pack( tpsit_d(jk,:,:)%vmp ,gsel )
    END DO
  ENDIF
!
! .. Sea ice 2D 
!
    DO jk=1,nt
        tzsit_r(jk,:)%esi = pack( tpsit(jk,:,:)%esi ,gsel )
        tzsit_r(jk,:)%asn = pack( tpsit(jk,:,:)%asn ,gsel )
        tzsit_r(jk,:)%fsi = pack( tpsit(jk,:,:)%fsi ,gsel )
        tzsit_r(jk,:)%hsi = pack( tpsit(jk,:,:)%hsi ,gsel )
        tzsit_r(jk,:)%hsn = pack( tpsit(jk,:,:)%hsn ,gsel )
        tzsit_r(jk,:)%rsn = pack( tpsit(jk,:,:)%rsn ,gsel )
        tzsit_r(jk,:)%tsf = pack( tpsit(jk,:,:)%tsf ,gsel )
        tzsit_r(jk,:)%ssi = pack( tpsit(jk,:,:)%ssi ,gsel )
        tzsit_r(jk,:)%age = pack( tpsit(jk,:,:)%age ,gsel )
        tzsit_r(jk,:)%vmp = pack( tpsit(jk,:,:)%vmp ,gsel )
    END DO
!
! .. Sea ice 3D
!
    DO jl=1,nl
      DO jk=1,nt
          tzsil_r(jl,jk,:)%ent = pack( tpsil(jl,jk,:,:)%ent ,gsel )
      END DO
    END DO 
! 
!
! 2.4. Update domain surface
! ---------------------------
!
    xdomsrf_r = sum( tzdom_r(:)%dxc*tzdom_r(:)%dyc )
!
!
! 2.5. Thermodynamics on the reduced grid
! ----------------------------------------
!
    IF ( nthermo==1 ) THEN
      CALL glt_thermo_r  &
        ( tzdom_r,zustar_r,tzmxl_r,tzatm_r,tzblkw_r,tzblki_r,tzbud_r,tzdia_r,  &
        tztfl_r,tzsit_r,tzsil_r )
    ENDIF
!
!
! 2.6. Apply the constraint
! --------------------------
!
    IF ( ntd==1 ) THEN
      CALL glt_constrain_r( tzdom_r,tzmxl_r,tzsit_r,tzsil_r,tzdia_r,tzsit_d_r )
    ENDIF
!
!
! 2.7. From reduced to global grid arrays
! ----------------------------------------
!
! .. Energy budget 
! Nothing to do : this is a non-cumulative diagnostic (initialized at 
! every glt_updtfl CALL)
!
! .. Diagnostics
! Note that only INTENT(inout) arguments must be unpacked
!
    tzdia0(:,:)%uvl = 0.
    tzdia0(:,:)%vvl = 0.
    tzdia0(:,:)%aiw = 0.
    tzdia0(:,:)%asi = 0.
    tzdia0(:,:)%amp = 0.
    tzdia0(:,:)%asn = 0.
    tzdia0(:,:)%cgl = 0.
    tzdia0(:,:)%dsa = 0.
    tzdia0(:,:)%dsn = 0.
    tzdia0(:,:)%dsi = 0.
    tzdia0(:,:)%dci = 0.
    tzdia0(:,:)%cst = 0.
    tzdia0(:,:)%dwi = 0.
    tzdia0(:,:)%lip = 0.
    tzdia0(:,:)%lsi = 0.
    tzdia0(:,:)%mrb = 0.
    tzdia0(:,:)%mrt = 0.
    tzdia0(:,:)%mrl = 0.
    tzdia0(:,:)%sie = 0.
    tzdia0(:,:)%sne = 0.
    tzdia0(:,:)%sni = 0.
    tzdia0(:,:)%snm = 0.
    tzdia0(:,:)%snml = 0.
    tzdia0(:,:)%sop = 0.
    tzdia0(:,:)%the = 0.
    tzdia0(:,:)%tin = 0.
    tzdia0(:,:)%tiw = 0.
    tzdia0(:,:)%qoi = 0.
    tzdia0(:,:)%xtr = 0.
    tzdia0(:,:)%ytr = 0.
    tzdia0(:,:)%sp1 = 0.
    tzdia0(:,:)%sp2 = 0.
    tzdia0(:,:)%sui = 0.
    tzdia0(:,:)%sut = 0.
    tzdia0(:,:)%sus = 0.
    tzdia0(:,:)%suw = 0.
    tzdia0(:,:)%sul = 0.
    tzdia0(:,:)%s_pr = 0.
    tzdia0(:,:)%s_prsn = 0.
    tzdia0(:,:)%o_pr = 0.
    tzdia0(:,:)%o_prsn = 0.
    tzdia0(:,:)%l_pr = 0.
    tzdia0(:,:)%l_prsn = 0.
    tzdia0(:,:)%subcio = 0.
    tzdia0(:,:)%snicio = 0.
    tzdia0(:,:)%hsicio = 0.
    tzdia0(:,:)%lmlcio = 0.
    tzdia0(:,:)%salcio = 0.
!
    tpdia%uvl = unpack( tzdia_r%uvl,gsel,tzdia0%uvl )
    tpdia%vvl = unpack( tzdia_r%vvl,gsel,tzdia0%vvl )
    tpdia%aiw = unpack( tzdia_r%aiw,gsel,tzdia0%aiw )
    tpdia%asi = unpack( tzdia_r%asi,gsel,tzdia0%asi )
    tpdia%amp = unpack( tzdia_r%amp,gsel,tzdia0%amp )
    tpdia%asn = unpack( tzdia_r%asn,gsel,tzdia0%asn )
    tpdia%cgl = unpack( tzdia_r%cgl,gsel,tzdia0%cgl )
    tpdia%dsa = unpack( tzdia_r%dsa,gsel,tzdia0%dsa )
    tpdia%dsn = unpack( tzdia_r%dsn,gsel,tzdia0%dsn )
    tpdia%dsi = unpack( tzdia_r%dsi,gsel,tzdia0%dsi )
    tpdia%dci = unpack( tzdia_r%dci,gsel,tzdia0%dci )
    tpdia%cst = unpack( tzdia_r%cst,gsel,tzdia0%cst )
    tpdia%dwi = unpack( tzdia_r%dwi,gsel,tzdia0%dwi )
    tpdia%lip = unpack( tzdia_r%lip,gsel,tzdia0%lip )
    tpdia%lsi = unpack( tzdia_r%lsi,gsel,tzdia0%lsi )
    tpdia%mrb = unpack( tzdia_r%mrb,gsel,tzdia0%mrb )
    tpdia%mrt = unpack( tzdia_r%mrt,gsel,tzdia0%mrt )
    tpdia%mrl = unpack( tzdia_r%mrl,gsel,tzdia0%mrl )
    tpdia%sie = unpack( tzdia_r%sie,gsel,tzdia0%sie )
    tpdia%sne = unpack( tzdia_r%sne,gsel,tzdia0%sne )
    tpdia%sni = unpack( tzdia_r%sni,gsel,tzdia0%sni )
    tpdia%snm = unpack( tzdia_r%snm,gsel,tzdia0%snm )
    tpdia%snml = unpack( tzdia_r%snml,gsel,tzdia0%snml )
    tpdia%sop = unpack( tzdia_r%sop,gsel,tzdia0%sop )
    tpdia%the = unpack( tzdia_r%the,gsel,tzdia0%the )
    tpdia%tin = unpack( tzdia_r%tin,gsel,tzdia0%tin )
    tpdia%tiw = unpack( tzdia_r%tiw,gsel,tzdia0%tiw )
    tpdia%qoi = unpack( tzdia_r%qoi,gsel,tzdia0%qoi )
    tpdia%xtr = unpack( tzdia_r%xtr,gsel,tzdia0%xtr )
    tpdia%ytr = unpack( tzdia_r%ytr,gsel,tzdia0%ytr )
    tpdia%sp1 = unpack( tzdia_r%sp1,gsel,tzdia0%sp1 )
    tpdia%sp2 = unpack( tzdia_r%sp2,gsel,tzdia0%sp2 )
    tpdia%sus = unpack( tzdia_r%sus,gsel,tzdia0%sus )
    tpdia%sut = unpack( tzdia_r%sut,gsel,tzdia0%sut )
    tpdia%sui = unpack( tzdia_r%sui,gsel,tzdia0%sui )
    tpdia%suw = unpack( tzdia_r%suw,gsel,tzdia0%suw )
    tpdia%s_pr = unpack( tzdia_r%s_pr,gsel,tzdia0%s_pr )
    tpdia%s_prsn = unpack( tzdia_r%s_prsn,gsel,tzdia0%s_prsn )
    tpdia%o_pr = unpack( tzdia_r%o_pr,gsel,tzdia0%o_pr )
    tpdia%o_prsn = unpack( tzdia_r%o_prsn,gsel,tzdia0%o_prsn )
    tpdia%subcio = unpack( tzdia_r%subcio,gsel,tzdia0%subcio )
    tpdia%snicio = unpack( tzdia_r%snicio,gsel,tzdia0%snicio )
    tpdia%hsicio = unpack( tzdia_r%hsicio,gsel,tzdia0%hsicio )
    tpdia%lmlcio = unpack( tzdia_r%lmlcio,gsel,tzdia0%lmlcio )
    tpdia%salcio = unpack( tzdia_r%salcio,gsel,tzdia0%salcio )
    tpdia%l_pr = unpack( tzdia_r%l_pr,gsel,tzdia0%l_pr )
    tpdia%l_prsn = unpack( tzdia_r%l_prsn,gsel,tzdia0%l_prsn )
    tpdia%sul = unpack( tzdia_r%sul,gsel,tzdia0%sul )
!
! .. Fluxes transmitted to the mixed layer
!
! Components of tptfl to be updated out of the reduced grid:
!  - llo : short wave on leads : only thermo_lead
!  - tlo : non solar on leads : only in thermo_lead
!  - wlo : fresh water mass flux through leads : only in thermo_lead
! (Redundant with what is in imod_thermo_lead... Caution !)
!
! Short wave on leads
    tztfl0(:,:)%llo = tptfl(:,:)%llo + tpblkw(:,:)%swa
    tptfl%llo = unpack( tztfl_r%llo,gsel,tztfl0%llo )
!
! Non solar flux on leads
    IF ( nsnwrad==1 ) THEN
        zsnflx(:,:) = -xmhofusn0*tpatm(:,:)%sop
      ELSE
        zsnflx(:,:) = 0.
    ENDIF
    tztfl0(:,:)%tlo = tptfl(:,:)%tlo +  &
      tpmxl(:,:)%qml + tpblkw(:,:)%nsf + zsnflx(:,:)
    tptfl%tlo = unpack( tztfl_r%tlo,gsel,tztfl0%tlo )
!
! Water flux on leads
    tztfl0(:,:)%wlo = tptfl(:,:)%wlo +  &
      tpatm(:,:)%lip + tpatm(:,:)%sop + tpblkw(:,:)%eva
    tptfl%wlo = unpack( tztfl_r%wlo,gsel,tztfl0%wlo )
!
! Short wave flux on leads
    tztfl0(:,:)%lio = 0.
    tptfl%lio = unpack( tztfl_r%lio,gsel,tztfl0%lio )
!
! Non solar flux under sea ice
    tztfl0(:,:)%tio = 0.
    tptfl%tio = unpack( tztfl_r%tio,gsel,tztfl0%tio )
!
! Salinity flux under sea ice
    tztfl0(:,:)%sio = 0.
    tptfl%sio = unpack( tztfl_r%sio,gsel,tztfl0%sio )
!
! Concentration/dilution flux under sea ice
    tztfl0(:,:)%cio = 0.
    tptfl%cio = unpack( tztfl_r%cio,gsel,tztfl0%cio )
!
! Water flux mass flux under sea ice
    tztfl0(:,:)%wio = 0.
    tptfl%wio = unpack( tztfl_r%wio,gsel,tztfl0%wio )
!
! Sea ice-ocean stress X-component
    tztfl0(:,:)%xio = 0.
    tptfl%xio = unpack( tztfl_r%xio,gsel,tztfl0%xio )
!
! Sea ice-ocean stress Y-component
    tztfl0(:,:)%yio = 0.
    tptfl%yio = unpack( tztfl_r%yio,gsel,tztfl0%yio )
!
! .. Sea ice 2D
!
    tzsit0(:,:,:)%esi = .false.
    tzsit0(:,:,:)%asn = albw
    tzsit0(:,:,:)%fsi = 0.
    tzsit0(:,:,:)%hsi = 0.
    tzsit0(:,:,:)%hsn = 0.
    tzsit0(:,:,:)%rsn = rhosnwmin
    tzsit0(:,:,:)%tsf = spread( tpmxl(:,:)%mlf,1,nt )
    tzsit0(:,:,:)%ssi = 0.
    tzsit0(:,:,:)%age = 0.
    tzsit0(:,:,:)%vmp = 0.
!
    DO jk=1,nt
      tpsit(jk,:,:)%esi = unpack( tzsit_r(jk,:)%esi,gsel,tzsit0(jk,:,:)%esi )
      tpsit(jk,:,:)%asn = unpack( tzsit_r(jk,:)%asn,gsel,tzsit0(jk,:,:)%asn )
      tpsit(jk,:,:)%fsi = unpack( tzsit_r(jk,:)%fsi,gsel,tzsit0(jk,:,:)%fsi )
      tpsit(jk,:,:)%hsi = unpack( tzsit_r(jk,:)%hsi,gsel,tzsit0(jk,:,:)%hsi )
      tpsit(jk,:,:)%hsn = unpack( tzsit_r(jk,:)%hsn,gsel,tzsit0(jk,:,:)%hsn )
      tpsit(jk,:,:)%rsn = unpack( tzsit_r(jk,:)%rsn,gsel,tzsit0(jk,:,:)%rsn )
      tpsit(jk,:,:)%tsf = unpack( tzsit_r(jk,:)%tsf,gsel,tzsit0(jk,:,:)%tsf )
      tpsit(jk,:,:)%ssi = unpack( tzsit_r(jk,:)%ssi,gsel,tzsit0(jk,:,:)%ssi )
      tpsit(jk,:,:)%age = unpack( tzsit_r(jk,:)%age,gsel,tzsit0(jk,:,:)%age )
      tpsit(jk,:,:)%vmp = unpack( tzsit_r(jk,:)%vmp,gsel,tzsit0(jk,:,:)%vmp )
    END DO
!
! .. Sea ice 3D
!
    tzsil0(:,:,:,:)%ent = 0.
!
    DO jl=1,nl
      DO jk=1,nt
        tpsil(jl,jk,:,:)%ent = unpack( tzsil_r(jl,jk,:)%ent,gsel,  &
          tzsil0(jl,jk,:,:)%ent )
      END DO
    END DO
!
!
!
! 3. Deallocations
! =================
!
! .. Note that on some platforms, in contradiction with the Fortran 95 norm,
! locally-allocated arrays are not released automatically. More over, some
! platforms require that the last allocated array in deallocated first.
! These rules have to be respected if more arrays are allocated/deallocated 
! here.
!
!    DEALLOCATE( ind_r )
    DEALLOCATE( tzsil_r )
    DEALLOCATE( tzsit_r )
    DEALLOCATE( tztfl_r )
    DEALLOCATE( tzdia_r )
    DEALLOCATE( tzbud_r )
    DEALLOCATE( tzblki_r )
    DEALLOCATE( tzblkw_r )
    DEALLOCATE( tzatm_r )
    DEALLOCATE( tzmxl_r )
    DEALLOCATE( zustar_r )
    DEALLOCATE( tzdom_r )
  ENDIF
!
!
!
! 4. Formats + Farewell message
! ==============================
!
1000 FORMAT( " Processor ", i5," ==> Running on ", i5,  &
  " points instead of ", i6, "(" i5, " times " , i5, ")" ) 
!
  IF (lp1) THEN
    WRITE(noutlu,*) ' '
    WRITE(noutlu,*) ' *** LEVEL 3 - END SUBROUTINE THERMO'
    WRITE(noutlu,*) ' '
  ENDIF
!
END SUBROUTINE glt_thermo

! ----------------------- END SUBROUTINE glt_thermo -------------------------
! -----------------------------------------------------------------------