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,tpglt, psst, psss, &
         pfsic, pfsit, psic, ptice, pice_alb)
!     #######################################################################
!
!!****  *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_types_glt , ONLY : t_glt
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)
TYPE(t_glt)         ,INTENT(INOUT)    :: tpglt     ! Gelato state variable
REAL, DIMENSION(:)  ,INTENT(IN)       :: psst      ! sea surface temperature (K)
REAL, DIMENSION(:)  ,INTENT(IN)       :: psss      ! sea surface salinity (psu)
REAL, DIMENSION(:)  ,INTENT(IN)       :: pfsic ! sea ice cover constraint ([0-1]). 
REAL, DIMENSION(:)  ,INTENT(IN)       :: pfsit ! sea ice thickness constraint (m). 
!
REAL, DIMENSION(:)  ,INTENT(OUT)      :: psic      ! Sea-ice Cover ([0-1])
REAL, DIMENSION(:)  ,INTENT(OUT)      :: ptice     ! Sea-ice temperature (K)
REAL, DIMENSION(:)  ,INTENT(OUT)      :: pice_alb  ! Sea-ce albedo ([0-1])
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(  SIZE(PSST)) :: zsst  ! sea surface temperature with frozen sea set at freezing point
REAL, DIMENSION(  SIZE(PSST)) :: zsic  ! Will hold a forcing SIC field, even if PFSIC 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(psss)
!
! 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 
tpglt%oce_all(:,1)%sml=psss(:)

! 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(psss,(/SIZE(psss),1/))) + xtt,(/SIZE(psss)/))

! Then replace freezing temp with Surfex-provided SST (PSST) 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=pfsic
   WHERE ( zsic(:) < 1.e-10 .AND. psst(:) > zsst(:) ) 
      zsst(:)=psst(:)
   ENDWHERE
ELSE 
   ! Implicit sea-ice cover
   WHERE (psst(:) - xtt > s%XFREEZING_SST + 0.1 )
      zsst(:)= psst(:) 
      zsic(:)=0.
   ELSEWHERE
      zsic(:)=1.
   ENDWHERE
ENDIF
!
!        Run Gelato for the length of the atmospheric time step
!________________________________________________________________________
!
!  Reset output accumulation/averaging fields 
psic    = 0.
ptice   = 0.
pice_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(psss) > 0) THEN 
      tpglt%oce_all(:,1)%tml=zsst(:)
      IF (s%LINTERPOL_SIC) tpglt%sit_d(1,:,1)%fsi=zsic(:)
      IF (s%LINTERPOL_SIT) tpglt%sit_d(1,:,1)%hsi=pfsit(:)
      ! Gelato will compute heat flux from ocean by itself, thanks to 
      ! imposed namelist parameter nextqoc=0 
      tpglt%oce_all(:,1)%qoc=0. 
      ! Zero Frazil flux
      tpglt%oce_all(:,1)%qml=0.
      ! Don't bother for sea level variations
      tpglt%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
      ! 
      tpglt%atm_all(:,1)%lip=s%XCPL_SEA_RAIN(:) / ptstep
      tpglt%atm_all(:,1)%sop=s%XCPL_SEA_SNOW(:) / ptstep
      ! Fluxes over Sea water
      tpglt%atm_wat(:,1)%eva=s%XCPL_SEA_EVAP(:) / ptstep
      tpglt%atm_wat(:,1)%swa=s%XCPL_SEA_SNET(:) / ptstep 
      tpglt%atm_wat(:,1)%nsf=s%XCPL_SEA_HEAT(:) / ptstep 
      tpglt%atm_wat(:,1)%dfl=s%XSI_FLX_DRV ! W m-2 K-1    
      ! Fluxes over Sea ice
      tpglt%atm_ice(1,:,1)%eva=s%XCPL_SEAICE_EVAP(:) / ptstep
      tpglt%atm_ice(1,:,1)%swa=s%XCPL_SEAICE_SNET(:) / ptstep 
      tpglt%atm_ice(1,:,1)%nsf=s%XCPL_SEAICE_HEAT(:) / ptstep 
      tpglt%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( tpglt%oce_all,tpglt%tml )
      CALL glt_getatmf( tpglt )
      CALL gltools_chkinp( 20010101,tpglt )
      !
      ! Compute gelato time index
      !
      tpglt%IND%CUR = ( ptimec + jt * xtstep ) / xtstep 
      !
      !       Let Gelato thermodynamic scheme run
      !
      CALL glt_gelato( tpglt )
      !
      ! Have Gelato feed its coupling ouptut interface
      ! 
      CALL glt_sndatmf( tpglt )
      CALL glt_sndmlrf( tpglt%bat,tpglt%dom,tpglt%atm_all,tpglt%tml, &
           tpglt%dia,tpglt%sit,tpglt%tfl,tpglt%ust,tpglt%all_oce )
      CALL wridia_ar5( tpglt )
      CALL gltools_chkout( 20010101,tpglt ) ! Does not actually work with Arpege
      ! Sum output fields over Gelato model time step duration
      psic     = psic     + tpglt%ice_atm(1,:,1)%fsi * xtstep
      ptice    = ptice    + tpglt%ice_atm(1,:,1)%tsf * xtstep
      pice_alb = pice_alb + tpglt%ice_atm(1,:,1)%alb * xtstep
   ENDIF
END DO
!   Average output fields over coupling time
psic     = psic     / (it * xtstep)
ptice    = ptice    / (it * xtstep) 
pice_alb = pice_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