seaice_gelato1dn.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.
!     #########
    SUBROUTINE seaice_gelato1d_n (S, HPROGRAM, PTIMEC, PTSTEP)
!     #######################################################################
!
!!****  *SEAICE_GELATO1D_n*  
!!
!!    PURPOSE
!!    -------
!     Run Gelato Sea-ice model, which provides sea-ice cover, temperature
!     and albedo
!     
!!**  METHOD
!!    ------
!     At relevant time steps :
!     i) Feed Gelato input interface structures with relevant variables 
!        taken in XCPL_xxx fields from MODD_SEAFLUX_n
!     ii) call the minimal Gelato sequence , namely : 
!        a) ingest and transform inputs regarding atm and sea state, 
!        b) run Gelato , 
!        c) have it process its ouptuts for atm and oce, 
!     iii) feed output arguments using content of Gelato output interface 
!        structure
!
!     Note : for the time being, SST and SSS are not averaged over coupling 
!     period
!
!!    EXTERNAL :
!!    ---------
!!    - a number of glt<xxx> routines
!!
!!    IMPLICIT ARGUMENTS :  
!!    ------------------- 
!!    Variables XCPL_xxx in MODD_SEAFLUX_n, and also LINTERPOL_SIC and LINTERPOL_SIT
!!      
!!    REFERENCE : 
!!    ---------
!!    Salas y Melia D (2002) A global coupled sea ice-ocean model. 
!!    Ocean Model 4:137-172
!!
!!    AUTHOR
!!    ------
!!     S.Senesi  *Meteo-France* CNRM - GAME
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     01/2014
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_seaflux_n, ONLY : seaflux_t
!
USE modd_csts,ONLY : xtt
USE modd_surf_par,   ONLY : xundef
USE modd_glt_param , ONLY : xtstep=>dtt, lwg, lp1, lp2, lp3, lp4, lp5, &
                            nprinto, gelato_dim=>nx

USE modi_glt_gelato
USE modi_glt_sndatmf
USE modi_glt_sndmlrf
USE modi_glt_getmlrf
USE modi_glt_getatmf
USE modi_gltools_chkinp
USE modi_gltools_chkout
USE mode_glt_stats
USE mode_gltools_swfrzt
USE mode_glt_dia_ar5
!
!
USE modi_get_luout
USE modi_abor1_sfx
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
TYPE(seaflux_t), INTENT(INOUT) :: S
!
 CHARACTER(LEN=6),    INTENT(IN)       :: HPROGRAM  ! program calling surf. schemes
REAL,                INTENT(IN)       :: PTIMEC    ! current duration since start of the run (s)
REAL,                INTENT(IN)       :: PTSTEP    ! surface time-step (s)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(  SIZE(S%XSST)) :: ZSST  ! sea surface temperature with frozen sea set at freezing point
REAL, DIMENSION(  SIZE(S%XSST)) :: ZSIC  ! Will hold a forcing SIC field, even if S%XFSIC is missing
!
INTEGER :: IT      ! total number of Gelato timesteps in one atmospheric timestep
INTEGER :: JT      ! Running index 1 -> IT
REAL    :: ZT      ! total number of Gelato timesteps in one atmospheric timestep
REAL    :: ZTIMEG  ! Gelato simulation time (s)

! The time step for ice coupling could be a namelist parameter
! for now, it will be set, here, as equal to the ice model time step
REAL :: ZICE_COUPLING_TSTEP
!
!
INTEGER         :: ILUOUT              ! output listing logical unit
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('SEAICE_GELATO1D',0,zhook_handle)
!
 CALL get_luout(hprogram,iluout)
!
! Must restore Gelato problem size (nx) to the correct value for the NPROMA block
! 
gelato_dim=SIZE(s%XSSS)
!
! Time steps stuff : default Gelato time step equals surface time step
!
IF (s%XSEAICE_TSTEP == xundef) THEN
  it=1
  zt=1.
  xtstep=ptstep
!
!* case of a Gelato time-step specified using NAM_SEAICE
ELSE
   IF ( ptstep < s%XSEAICE_TSTEP ) THEN
      CALL abor1_sfx("XSEAICE_TSTEP SHOULD BE EQUAL OR LESS THAN ATMOSPHERIC TIME STEP")
   ELSE
      it=nint(ptstep/s%XSEAICE_TSTEP)
      zt=float(it)
      xtstep=ptstep/zt
   ENDIF
ENDIF
!
!          Initializations
!________________________________________________________________________
!
! Ocean part 
!----------------------------------------------------------------------------------
! Surface salinity 
s%TGLT%oce_all(:,1)%sml=s%XSSS(:)

! Ensure that SSS-dependant freezing point temperature is used on
! locations where SIC (or SST) forcing value calls for it

! First init ZSST with freezing temperature (which depends on local salinity)
! (Gelato uses Celsius scale for freezing temperatures)
zsst=reshape(glt_swfrzt2d(reshape(s%XSSS,(/SIZE(s%XSSS),1/))) + xtt,(/SIZE(s%XSSS)/))

! Then replace freezing temp with Surfex-provided SST (S%XSST) where
! there is no (explicit or implicit) seaice and temperature is warmer
! than freezing point. And inits ZSIC accordingly
IF (s%LINTERPOL_SIC) THEN  
   zsic=s%XFSIC
   WHERE ( zsic(:) < 1.e-10 .AND. s%XSST(:) > zsst(:) ) 
      zsst(:)=s%XSST(:)
   ENDWHERE
ELSE 
   ! Implicit sea-ice cover
   WHERE (s%XSST(:) - xtt > s%XFREEZING_SST + 0.1 )
      zsst(:)= s%XSST(:) 
      zsic(:)=0.
   ELSEWHERE
      zsic(:)=1.
   ENDWHERE
ENDIF
!
!        Run Gelato for the length of the atmospheric time step
!________________________________________________________________________
!
!  Reset output accumulation/averaging fields 
s%XSIC    = 0.
s%XTICE   = 0.
s%XICE_ALB= 0.
lp1 = (lwg.AND.nprinto>=1)
lp2 = (lwg.AND.nprinto>=2)
lp3 = (lwg.AND.nprinto>=3)
lp4 = (lwg.AND.nprinto>=4)
lp5 = (lwg.AND.nprinto>=5)

DO jt=1,it
   IF (SIZE(s%XSSS) > 0) THEN 
      s%TGLT%oce_all(:,1)%tml=zsst(:)
      IF (s%LINTERPOL_SIC) s%TGLT%sit_d(1,:,1)%fsi=zsic(:)
      IF (s%LINTERPOL_SIT) s%TGLT%sit_d(1,:,1)%hsi=s%XFSIT(:)
      ! Gelato will compute heat flux from ocean by itself, thanks to 
      ! imposed namelist parameter nextqoc=0 
      s%TGLT%oce_all(:,1)%qoc=0. 
      ! Zero Frazil flux
      s%TGLT%oce_all(:,1)%qml=0.
      ! Don't bother for sea level variations
      s%TGLT%oce_all(:,1)%ssh=0.
      ! (velocity components are useless in Surfex 1D setting)
      !
      ! Atmosphere part 
      !----------------
      ! Feed Gelato input structure with flux values from XCPL_xx
      ! 
      s%TGLT%atm_all(:,1)%lip=s%XCPL_SEA_RAIN(:) / ptstep
      s%TGLT%atm_all(:,1)%sop=s%XCPL_SEA_SNOW(:) / ptstep
      ! Fluxes over Sea water
      s%TGLT%atm_wat(:,1)%eva=s%XCPL_SEA_EVAP(:) / ptstep
      s%TGLT%atm_wat(:,1)%swa=s%XCPL_SEA_SNET(:) / ptstep 
      s%TGLT%atm_wat(:,1)%nsf=s%XCPL_SEA_HEAT(:) / ptstep 
      s%TGLT%atm_wat(:,1)%dfl=s%XSI_FLX_DRV ! W m-2 K-1    
      ! Fluxes over Sea ice
      s%TGLT%atm_ice(1,:,1)%eva=s%XCPL_SEAICE_EVAP(:) / ptstep
      s%TGLT%atm_ice(1,:,1)%swa=s%XCPL_SEAICE_SNET(:) / ptstep 
      s%TGLT%atm_ice(1,:,1)%nsf=s%XCPL_SEAICE_HEAT(:) / ptstep 
      s%TGLT%atm_ice(1,:,1)%dfl=s%XSI_FLX_DRV ! W m-2 K-1   
      ! (stress components are useless in Surfex 1D setting)
      !
      !       Let Gelato process its input data
      !
      CALL glt_getmlrf( s%TGLT%oce_all,s%TGLT%tml )
      CALL glt_getatmf( s%TGLT )
      CALL gltools_chkinp( 20010101,s%TGLT )
      !
      ! Compute gelato time index
      !
      s%TGLT%IND%CUR = ( ptimec + jt * xtstep ) / xtstep 
      !
      !       Let Gelato thermodynamic scheme run
      !
      CALL glt_gelato( s%TGLT )
      !
      ! Have Gelato feed its coupling ouptut interface
      ! 
      CALL glt_sndatmf( s%TGLT )
      CALL glt_sndmlrf( s%TGLT%bat,s%TGLT%dom,s%TGLT%atm_all,s%TGLT%tml, &
                        s%TGLT%dia,s%TGLT%sit,s%TGLT%tfl,s%TGLT%ust,s%TGLT%all_oce )
      CALL wridia_ar5( s%TGLT )
      CALL gltools_chkout( 20010101,s%TGLT ) ! Does not actually work with Arpege
      ! Sum output fields over Gelato model time step duration
      s%XSIC     = s%XSIC     + s%TGLT%ice_atm(1,:,1)%fsi * xtstep
      s%XTICE    = s%XTICE    + s%TGLT%ice_atm(1,:,1)%tsf * xtstep
      s%XICE_ALB = s%XICE_ALB + s%TGLT%ice_atm(1,:,1)%alb * xtstep
   ENDIF
END DO
!   Average output fields over coupling time
s%XSIC     = s%XSIC     / (it * xtstep)
s%XTICE    = s%XTICE    / (it * xtstep) 
s%XICE_ALB = s%XICE_ALB / (it * xtstep) 
!
! Resets input accumulation fields for next step 
!
s%XCPL_SEA_RAIN=0.
s%XCPL_SEA_SNOW=0.
s%XCPL_SEA_EVAP=0.
s%XCPL_SEA_SNET=0.
s%XCPL_SEA_HEAT=0.
s%XCPL_SEAICE_EVAP=0.
s%XCPL_SEAICE_SNET=0.
s%XCPL_SEAICE_HEAT=0.
!
IF (lhook) CALL dr_hook('SEAICE_GELATO1D',1,zhook_handle)
!!-------------------------------------------------------------------------------
!!-----------------------------------------------------------------------------
END SUBROUTINE seaice_gelato1d_n