modi_glt_updbud.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_updbud =======================
! =======================================================================
!
! Goal:
! -----
!   This module contains a subroutine that manages the energy budget 
! variable on the whole grid:
!       - total sea ice gltools_enthalpy (reference temperature is the melting
!       temperature of snow and ice)
!       - total sea ice latent heat (i.e. the latent heat needed to 
!       melt sea ice and snow covers on the whole ocean domain). Note 
!       sea ice stored heat is included in that.
!       - total stored heat in sea ice.
!       - input and glt_output heat concerning leads and sea ice.
!   All these quantities are given in J. The energy budget is also 
! printed (total values on the whole mesh).
!   Note that this routine tends to mix pure prints (can be disabled), but
! that part of the calculations are also used by the physics of Gelato.
! Before clean splitting of these functionalities is done, be cautious if
! you wish to modify the code !
!   Note also that some calculations are global, i.e. when glt_avg routine
! is used, the information from all processors must be gathered. Hence 
! these instructions should not be conditional to lp1 = .TRUE., but 
! rather to nprinto = 1 !
!
! Created : 2001/08 (D. Salas y Melia)
! Modified: 2009/12 (D. Salas y Melia) Enthalpy replaces temperature
!   as a state variable
! Modified: 2010/01 (D. Salas y Melia) Introduce fresh water budget
! Modified: 2010/03 (D. Salas y Melia) Introduce salinity budget
!
! --------------------- BEGIN MODULE modi_glt_updbud ------------------------
!
!THXS_SFX!MODULE modi_glt_updbud
!THXS_SFX!INTERFACE 
!THXS_SFX!!
!THXS_SFX!SUBROUTINE glt_updbud  &
!THXS_SFX!  ( kinit,omsg,tpdom,tpmxl,tptfl,tpatm,tpblkw,tpblki,tpsit,tpsil,tpbud )
!THXS_SFX!  USE modd_types_glt
!THXS_SFX!  USE modd_glt_param
!THXS_SFX!  INTEGER, INTENT(in) ::  &
!THXS_SFX!         kinit
!THXS_SFX!  CHARACTER(*), INTENT(in) ::  &
!THXS_SFX!        omsg
!THXS_SFX!  TYPE(t_dom), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpdom
!THXS_SFX!  TYPE(t_mxl), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpmxl
!THXS_SFX!  TYPE(t_tfl), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tptfl
!THXS_SFX!  TYPE(t_atm), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpatm
!THXS_SFX!  TYPE(t_blk), DIMENSION(nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpblkw
!THXS_SFX!  TYPE(t_blk), DIMENSION(nt,nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpblki
!THXS_SFX!  TYPE(t_sit), DIMENSION(nt,nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpsit
!THXS_SFX!  TYPE(t_vtp), DIMENSION(nl,nt,nx,ny), INTENT(in) ::  &
!THXS_SFX!        tpsil
!THXS_SFX!  TYPE(t_bud), DIMENSION(nx,ny), INTENT(inout) ::  &
!THXS_SFX!        tpbud
!THXS_SFX!END SUBROUTINE glt_updbud
!THXS_SFX!!
!THXS_SFX!END INTERFACE
!THXS_SFX!END MODULE modi_glt_updbud
!
! ----------------------- END MODULE modi_glt_updbud ------------------------
!
!
!
! -----------------------------------------------------------------------
! ------------------------- SUBROUTINE glt_updbud ---------------------------
!
! .. Subroutine used to check global heat budget.
!
SUBROUTINE glt_updbud  &
  ( kinit,omsg,tpdom,tpmxl,tptfl,tpatm,tpblkw,tpblki,tpsit,tpsil,tpbud )
!
  USE modd_types_glt
  USE modd_glt_param
  USE modd_glt_const_thm
  USE mode_glt_stats
  USE mode_glt_info
!
  IMPLICIT NONE
!
  INTEGER, INTENT(in) ::  &
         kinit
  CHARACTER(*), INTENT(in) ::  &
         omsg
  TYPE(t_dom), DIMENSION(nx,ny), INTENT(in) ::  &
        tpdom
  TYPE(t_mxl), DIMENSION(nx,ny), INTENT(in) ::  &
        tpmxl
  TYPE(t_tfl), DIMENSION(nx,ny), INTENT(in) ::  &
        tptfl
  TYPE(t_atm), DIMENSION(nx,ny), INTENT(in) ::  &
        tpatm
  TYPE(t_blk), DIMENSION(nx,ny), INTENT(in) ::  &
        tpblkw
  TYPE(t_blk), DIMENSION(nt,nx,ny), INTENT(in) ::  &
        tpblki
  TYPE(t_sit), DIMENSION(nt,nx,ny), INTENT(in) ::  &
        tpsit
  TYPE(t_vtp), DIMENSION(nl,nt,nx,ny), INTENT(in) ::  &
        tpsil
  TYPE(t_bud), DIMENSION(nx,ny), INTENT(inout) ::  &
        tpbud
!
  INTEGER ::  &
        jl
  REAL ::  &
        zenthalpy0,zhiit0,zbiit0,zhii0,zhio0,   &
        zhli0,zhlo0,zwio0,zwlo0,zwii0,zwli0,zcio0,zclo0,  &
        zwater0,zsalt0
  REAL, DIMENSION(nx,ny) ::  &
        zfsit
  REAL, DIMENSION(nx,ny) ::  &
        zenthalpy2,zhiit2,zbiit2,  &
        znli2,zniit2,zhli2,zhlo2,zwii2,zwli2,zwater2,zsalt2
  REAL, DIMENSION(nt,nx,ny) ::  &
        zenthalpys,zenthalpyi,zmsi,zmsn
!
!
!
! 1. Initializations
! ==================
!
  IF (lp1) THEN
      WRITE(noutlu,*) ' ' 
      WRITE(noutlu,*) ' **** glt_updbud ****' 
      WRITE(noutlu,*) omsg, '   (energy fluxes in W.m-2)'
      WRITE(noutlu,*) omsg, '   (water and salt fluxes in kg.m-2.day-1)'
  ENDIF
!
!
! 1.1. Initialize arrays
! -----------------------
!
  zfsit(:,:) = sum( tpsit(:,:,:)%fsi,dim=1 )
  zmsi(:,:,:) = rhoice * tpsit(:,:,:)%fsi * tpsit(:,:,:)%hsi
  zmsn(:,:,:) = tpsit(:,:,:)%rsn * tpsit(:,:,:)%fsi * tpsit(:,:,:)%hsn
!
!
! 1.2. Print information 
! -----------------------
!
  CALL glt_info_si( omsg,tpdom,tpsit=tpsit )
!
!
!
! 2. Incoming fluxes
! ===================
!
! 2.1. Incoming fluxes on sea ice (top+bottom of the slab)
! ---------------------------------------------------------
!
! .. Note that the incoming energy at the top represents on the one
! hand the sum of solar and non solar fluxes, the snow layer latent
! heat change due to new snowfalls, and the related snow layer gltools_enthalpy 
! variation.
! (note that tpatm%sop is in kg.m-2.s-1)
!
  IF ( kinit==1 ) THEN
!
! Energy
      zniit2(:,:) = tpatm(:,:)%sop *  &
         sum( tpsit(:,:,:)%fsi*tpsil(nl,:,:,:)%ent, dim=1 ) 
!        ( -xmhofusn0*zfsit(:,:) +  &
!          SUM( tpsit(:,:,:)%fsi *  &
!               ( cpice0*tpsit(:,:,:)%tsf ),DIM=1 ) )  ! +  &
!                    tpatm(:,:)%lip*  &
!          xmhofusn0*zfsit(:,:)
! atm -> ice energy
      zhiit2(:,:) =  &
        sum( tpsit(:,:,:)%fsi*  &
        ( tpblki(:,:,:)%nsf+tpblki(:,:,:)%swa ), dim=1 )
! oce -> ice energy ( qoc can only be determined a posteriori, in glt_updhsi_r )
      zbiit2(:,:) =  &
        tpmxl(:,:)%qml*zfsit(:,:)    ! + tpmxl(:,:)%qoc
!
      tpbud(:,:)%nii = zniit2(:,:)
      tpbud(:,:)%hii = zhiit2(:,:)+zniit2(:,:)
      tpbud(:,:)%bii = zbiit2(:,:)
!
! Water
      zwii2(:,:) = zfsit(:,:)* ( tpatm(:,:)%sop + tpatm(:,:)%lip ) +  &
        sum( tpsit(:,:,:)%fsi*tpblki(:,:,:)%eva, dim=1 )
      tpbud(:,:)%wii = zwii2(:,:)*xday2sec
  ENDIF
!
  IF ( nprinto>=1 ) THEN
!
! Energy
      zhiit0 = glt_avg( tpdom,tpbud(:,:)%hii,0 )
      zbiit0 = glt_avg( tpdom,tpbud(:,:)%bii,0 )
      zhii0 = zhiit0 + zbiit0
!
! Water
      zwii0 = glt_avg( tpdom,tpbud(:,:)%wii,0 )
!
      IF (lwg) THEN
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
        WRITE(noutlu,*) '    Incoming ENERGY under sea ice      :',  &
          zbiit0
        WRITE(noutlu,*) '    Incoming ENERGY top of sea ice     :',  &
          zhiit0
        WRITE(noutlu,*) '    Total incoming ENERGY on sea ice   :',  &
          zhii0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
        WRITE(noutlu,*) '    Total incoming WATER on sea ice    :',  &
          zwii0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
  ENDIF
!
!
! 2.2. Incoming fluxes on leads
! ------------------------------
!
  IF ( kinit==1 ) THEN
!
! Energy
      zhli2(:,:) =  &
        ( 1.-zfsit(:,:) )*   &
        ( tpmxl(:,:)%qml+tpblkw(:,:)%nsf+tpblkw(:,:)%swa )    
      IF ( nsnwrad==1 ) THEN
          znli2(:,:) = ( 1.-zfsit(:,:) )*  &
!            ( cpice0*tpmxl(:,:)%mlf-xmhofusn0 )*tpatm(:,:)%sop
            ( -xmhofusn0 )*tpatm(:,:)%sop
        ELSE
          znli2(:,:) = 0.
      ENDIF
      tpbud(:,:)%nli = znli2(:,:)
      tpbud(:,:)%hli = zhli2(:,:)+znli2(:,:)
!
! Water
      zwli2(:,:) = ( 1.-zfsit(:,:) )*  &
        ( tpatm(:,:)%sop + tpatm(:,:)%lip + tpblkw(:,:)%eva )*xday2sec
      tpbud(:,:)%wli = zwli2(:,:)
  ENDIF
!
  IF ( nprinto>=1 ) THEN
!
! Energy
      zhli0 = glt_avg( tpdom,tpbud(:,:)%hli,0 )
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total incoming ENERGY (on leads)   :',  &
          zhli0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
!
! Water
      zwli0 = glt_avg( tpdom,tpbud(:,:)%wli,0 )
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total incoming WATER (on leads)    :',  &
          zwli0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF 
  ENDIF
!
!
!
! 3. Outgoing fluxes
! ===================
!
! 3.1. Outgoing fluxes at the bottom of the sea ice slab
! -------------------------------------------------------
!
! Energy
  tpbud(:,:)%hio = tptfl(:,:)%lio+tptfl(:,:)%tio
!
  IF ( nprinto>=1 ) THEN
      zhio0 = glt_avg( tpdom,tpbud(:,:)%hio,0 )
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total outgoing ENERGY (under sea ice) :',  &
          zhio0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF 
!
! Water
! This flux is just described by tptfl%wio
!
      zwio0 = glt_avg( tpdom,tptfl(:,:)%wio,0 )*xday2sec
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total outgoing WATER (under sea ice)  :',  &
          zwio0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
!
! Salt
! A conversion is necessary to convert the concentration/dilution flux to
! salt mass changes ( in kg.m-2 ):
!     D S = sml * cio
!
      zsalt2(:,:) = tptfl(:,:)%sio
      zcio0 = glt_avg( tpdom,zsalt2,0 )*xday2sec
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total outgoing SALT (under sea ice)   :',  &
          zcio0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
  ENDIF
!
!
! 3.2. Outgoing fluxes through the water surface
! -----------------------------------------------
!
! Energy
  zhlo2(:,:) = tptfl(:,:)%llo+tptfl(:,:)%tlo 
!
  tpbud(:,:)%hlo = zhlo2(:,:)
!
  IF ( nprinto>=1 ) THEN
      zhlo0 = glt_avg( tpdom,tpbud(:,:)%hlo,0 )
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total outgoing ENERGY (under leads)   :',  &
          zhlo0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
!
! Water
      zwlo0 = glt_avg( tpdom,tptfl%wlo,0 )*xday2sec
      IF (lwg) THEN
        WRITE(noutlu,*) '    Total outgoing WATER (under leads)    :',  &
          zwlo0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
!
! Salt
! A conversion is necessary to convert the concentration/dilution flux to
! salt mass changes ( in kg.m-2 ):
!     D S = sml * wlo
!
!      zsalt2(:,:) = 1.e-3 * tpmxl(:,:)%sml * tptfl(:,:)%wlo
!      zclo0 = glt_avg( tpdom,zsalt2,0 )*xday2sec
!      WRITE(noutlu,*) '    Total outgoing SALT (under leads)     :',  &
!        zclo0
!      WRITE(noutlu,*)  &
!        '--------------------------------------------------------------------'
  ENDIF
!
!
!
! 4. Energy/water/salt retained by sea ice/snow per square meter
! ===============================================================
!
! 4.1. Enthalpy in sea ice/snow per square meter
! ----------------------------------------------
!
  zenthalpyi(:,:,:) = 0.
  DO jl=1,nilay
    zenthalpyi(:,:,:) = zenthalpyi(:,:,:) +  &
      sf3tinv(jl) * zmsi(:,:,:) * tpsil(jl,:,:,:)%ent
  END DO
!
  zenthalpys(:,:,:) =  &
    zmsn(:,:,:) * sum( tpsil(nilay+1:nl,:,:,:)%ent,dim=1 )/float(nslay)
!
  zenthalpy2(:,:) = sum( zenthalpys(:,:,:)+zenthalpyi(:,:,:), dim=1 )
!
! .. Initial gltools_enthalpy is set to computed gltools_enthalpy if kinit flag is on. 
! Else initial gltools_enthalpy is kept as such.
!
  IF ( kinit==1 ) THEN
      tpbud(:,:)%eni = zenthalpy2(:,:)
      tpbud(:,:)%enn = zenthalpy2(:,:)
    ELSE
      tpbud(:,:)%enn = zenthalpy2(:,:)
  ENDIF
!
  IF ( nprinto>=1 ) THEN
      zenthalpy0 = glt_avg( tpdom, tpbud(:,:)%enn-tpbud(:,:)%eni,0 ) / dtt
      IF (lwg) THEN
        WRITE(noutlu,*)  &
          '    D(enthalpy) from beg. of time step :',  &
          zenthalpy0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
  ENDIF
!
!
! 4.2. Water stored by sea ice and snow
! --------------------------------------
!
! .. Note that we exclude salt (only fresh water is taken into account)  
!
! Sea ice + snow fresh water content
  zwater2(:,:) =  xday2sec*( &
    sum( zmsi(:,:,:)*( 1.-1.e-3*tpsit(:,:,:)%ssi ), dim=1 ) +  &
    sum( zmsn(:,:,:), dim=1 ) )
!
! .. Initial fresh water content of sea ice + snow cover is initialised if
! kinit flag is on.
!
  IF ( kinit==1 ) THEN
      tpbud(:,:)%fwi = zwater2(:,:)
  ENDIF
  tpbud(:,:)%fwn = zwater2(:,:)
!
  IF ( nprinto>=1 ) THEN
      zwater0 = glt_avg( tpdom,tpbud(:,:)%fwn-tpbud(:,:)%fwi,0 ) / dtt
      IF (lwg) THEN
        WRITE(noutlu,*) '    D(water) from beg. of time step   :',  &
          zwater0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
  ENDIF
! 
!
! 4.3. Salt stored by sea ice
! ----------------------------
!
! Sea ice salt content
  zsalt2(:,:) = xday2sec*  &
    sum( zmsi(:,:,:)*1.e-3*tpsit(:,:,:)%ssi, dim=1 )
!
! .. Initial salt content of sea ice is initialised if kinit flag is on.
!
  IF ( kinit==1 ) THEN
      tpbud(:,:)%isi = zsalt2(:,:)
  ENDIF
  tpbud(:,:)%isn = zsalt2(:,:)
!
  IF ( nprinto>=1 ) THEN
      zsalt0 = glt_avg( tpdom,tpbud(:,:)%isn-tpbud(:,:)%isi,0 ) / dtt
      IF (lwg) THEN
        WRITE(noutlu,*) '    D(salt) from beg. of time step    :',  &
          zsalt0
        WRITE(noutlu,*)  &
          '--------------------------------------------------------------------'
      ENDIF
  ENDIF
!
!
! 5. Final prints: general budget
! ===============================
!
! .. At the end of the time step, the sum of incoming energy on leads
! and on sea ice should be equal to the sum of the total variation of
! gltools_enthalpy in sea ice+snow and of the total outgoing energy.
!    A print of these two quantities is done here.
!
  IF (lp1) THEN
      WRITE(noutlu,*) '    Total incoming energy (leads+ice)  :',  &
        zhii0+zhli0
      WRITE(noutlu,*) '    D(enthalpy) + outg. ENERGY         :',  &
        zenthalpy0+zhio0+zhlo0
      WRITE(noutlu,*) '    ENERGY BALANCE                     :',  &
        zenthalpy0+zhio0+zhlo0-(zhii0+zhli0)
      WRITE(noutlu,*)  &
        '--------------------------------------------------------------------'
      WRITE(noutlu,*) '    Total incoming + outgoing water (leads+ice)   :',  &
        zwii0+zwli0-zwio0-zwlo0
      WRITE(noutlu,*) '    WATER  BALANCE                                :',  &
        zwii0+zwli0-zwio0-zwlo0-zwater0
      WRITE(noutlu,*)  &
        '--------------------------------------------------------------------'
      WRITE(noutlu,*) '    SALT  BALANCE                                 :',  &
        -zcio0-zsalt0
      WRITE(noutlu,*)  &
        '--------------------------------------------------------------------'
  ENDIF
!
END SUBROUTINE glt_updbud

! ----------------------- END SUBROUTINE glt_updbud -------------------------
! -----------------------------------------------------------------------