read_seaicen.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 read_seaice_n (&
                                 sg, s, &
                                hprogram,klu,kluout)
!     #########################################
!
!!****  *READ_SEAICE_n* - read seaice scheme variables
!!
!!
!!    PURPOSE : feed seaice scheme state variable and 'domain' structure
!!    -------
!!
!!**  METHOD : 
!!    -------
!!      For now, only Gelato model is handled
!!
!!      for state variable : quite standard in Surfex : use READ_SURF with 
!!         relevant field names (same names as in genuine gelato restarts)
!!      for domain information : copy from MODD_SEAFLUX_GRID
!!      for bathymetry : copy from MODD_SEAFLUX
!!
!!    EXTERNALS : READ_SURF, GLTOOLS_ALLOC, GET_TYPE_DIM, ABOR1_SFX
!!    --------
!!
!!    IMPLICIT ARGUMENTS : Gelato state variable, and a few namelist parameters
!!    ------------------
!!
!!    REFERENCE : routine restartr in original Gelato sources (V6.0.20)
!!    ---------
!!
!!    AUTHOR : S. Sénési   *Meteo France*       
!!    ------
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2014
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!
!
!
!
USE modd_seaflux_grid_n, ONLY : seaflux_grid_t
USE modd_seaflux_n, ONLY : seaflux_t
!
USE modd_csts, ONLY           : xpi, xttsi, xtt
USE modd_surf_par,       ONLY : xundef
USE modd_sfx_oasis,      ONLY : lcpl_seaice
USE modd_water_par,      ONLY : xalbseaice
!
USE modd_types_glt,   ONLY : t_glt
USE modd_glt_param, ONLY : nl, nt, nx, ny, nxglo, nyglo, xdomsrf, &
                           xdomsrf_g, nprinto, cfsidmp, chsidmp,  &
                           xfsidmpeft, xhsidmpeft, ntd
USE modd_glt_const_thm, ONLY : epsil1
USE lib_mpp,            ONLY : mpp_sum
USE modi_glt_sndatmf
USE modi_gltools_alloc
USE modi_gltools_readnam
!
USE modi_read_surf
USE modi_interpol_sst_mth
!
USE modi_get_luout
USE modi_abor1_sfx
USE modd_surfex_omp, ONLY : nblocktot
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
!
!
TYPE(seaflux_grid_t), INTENT(INOUT) :: sg
TYPE(seaflux_t), INTENT(INOUT) :: s
!
 CHARACTER(LEN=6),  INTENT(IN)  :: hprogram ! calling program
INTEGER,           INTENT(IN)  :: klu      ! number of sea patch point
INTEGER,           INTENT(IN)  :: kluout
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER           :: iresp          ! Error code after reading
!
INTEGER           :: jmth, inmth
 CHARACTER(LEN=2 ) :: ymth
 CHARACTER(LEN=5)  :: ylvl
!
 CHARACTER(LEN=12) :: ycateg         ! category to read
 CHARACTER(LEN=12) :: ylevel         ! Level to read
 CHARACTER(LEN=12) :: yrecfm         ! Name of the article to be read
 CHARACTER(LEN=200) :: ymess         ! Error Message
!
INTEGER :: jx,jk,jl                 ! loop counter on ice categories and layers and grid points
INTEGER :: inl_in_file,int_in_file  ! file values for ice catgories and layers numbers
REAL :: zfsit
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('READ_SEAICE_n',0,zhook_handle)
!
IF (.NOT.s%LHANDLE_SIC) THEN 
   ALLOCATE(s%XSIC(0))
   IF (lhook) CALL dr_hook('READ_SEAICE_n',1,zhook_handle)
   RETURN
ENDIF
!
nx=klu
!
ALLOCATE(s%XSIC(klu))
s%XSIC(:)=xundef
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Dealing with external sea-ice cover data (either for nudging or forcing)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IF(s%LINTERPOL_SIC)THEN
   !
   ALLOCATE(s%XFSIC(klu))
   !
   !Precedent, Current, Next, and Second-next Monthly SIC
   inmth=4   
   !
   ALLOCATE(s%XSIC_MTH(klu,inmth))
   DO jmth=1,inmth
      WRITE(ymth,'(I2)') (jmth-1)
      yrecfm='SIC_MTH'//adjustl(ymth(:len_trim(ymth)))
      CALL read_surf(&
                     hprogram,yrecfm,s%XSIC_MTH(:,jmth),iresp)
      CALL check_seaice(yrecfm,s%XSIC_MTH(:,jmth))
   ENDDO
   !
   CALL interpol_sst_mth(s, &
                         s%TTIME%TDATE%YEAR,s%TTIME%TDATE%MONTH,s%TTIME%TDATE%DAY,'C',s%XFSIC)
   !
   IF (any(s%XFSIC(:)>1.0).OR.any(s%XFSIC(:)<0.0)) THEN
     CALL abor1_sfx('READ_SEAICE_n: FSIC should be >=0 and <=1') 
   ENDIF                 
   !
ELSE
   ! 
   ALLOCATE(s%XFSIC(0))
   ALLOCATE(s%XSIC_MTH(0,0))
   !
ENDIF
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!                             Worrying about a seaice scheme 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
 CALL read_surf(&
                     hprogram,'SEAICE_SCHEM',s%CSEAICE_SCHEME,iresp)
!
IF (trim(s%CSEAICE_SCHEME) == 'GELATO' .AND. (nblocktot>1)) THEN 
    CALL abor1_sfx("READ_SEAICE_n: GELATO CANNOT YET RUN MULTI-THREAD")
ENDIF
!
IF (trim(s%CSEAICE_SCHEME) == 'NONE' ) THEN
   IF (s%LINTERPOL_SIC ) THEN
      s%XTICE=s%XSST
      s%XSIC=s%XFSIC
      s%XICE_ALB=xalbseaice           
      IF (lhook) CALL dr_hook('READ_SEAICE_n',1,zhook_handle)
      RETURN
   ELSE
      CALL abor1_sfx("READ_SEAICE_n: MUST HAVE CINTERPOL_SIC /= NONE WITH CSEAICE_SCHEME == NONE ") 
   ENDIF
ELSE 
   IF (trim(s%CSEAICE_SCHEME) /= 'GELATO') THEN 
      WRITE(kluout,*)'READ_SEAICE_n:CSEAICE_SCHEME read in PREP, ',s%CSEAICE_SCHEME,', is not yet handled'
      CALL abor1_sfx("READ_SEAICE_n:CAN ONLY HANDLE GELATO SEAICE MODEL YET (and not the one quoted in PREP)") 
   ENDIF
ENDIF
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!                             Start of Gelato specifics
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IF(lcpl_seaice)THEN       
   CALL abor1_sfx('READ_SEAICEN: CANNOT YET MANAGE BOTH TRUE LCPL_SEAICE AND CSEAICE_SCHEME = GELATO')
ENDIF
nxglo=nx
#if ! defined in_arpege
 CALL mpp_sum(nxglo) ! Should also sum up over NPROMA blocks, in Arpege; but not that easy....
#else
IF (nprinto > 0) THEN
   WRITE(kluout,*)'Gelato cannot yet compute global averages when running in Arpege (because of collective comm vs. NPROMA blocks)'
ENDIF
nxglo=max(nxglo,1)
#endif
!
! Use convention XSIC_EFOLDING_TIME=0 for avoiding any relaxation 
! toward SIC observation and impose it.
!
IF(s%LINTERPOL_SIC)THEN
  IF (s%XSIC_EFOLDING_TIME==0.0) THEN 
     cfsidmp='PRESCRIBE'
  ELSE
     cfsidmp='DAMP'
     xfsidmpeft=s%XSIC_EFOLDING_TIME 
  ENDIF
ENDIF
!
IF(s%LINTERPOL_SIT)THEN
  IF (s%XSIT_EFOLDING_TIME==0.0) THEN 
     chsidmp='PRESCRIBE'
  ELSE
     chsidmp='DAMP_FAC'
     xhsidmpeft= s%XSIT_EFOLDING_TIME 
  ENDIF
ENDIF
!
!* Physical dimensions are set for Gelato , as a 1D field (second dimension is degenerated)
!
! Supersedes Gelato hard defaults with a Gelato genuine namelist 
! if available (for Gelato wizzards !)
 CALL gltools_readnam(.false.,kluout)  
!
ny=1
nyglo=1
 CALL gltools_alloc(s%TGLT)
!
!*       0.     Check dimensions : number of layers and ice categories
!
 CALL read_surf(&
                     hprogram,'ICENL',inl_in_file,iresp)
IF (inl_in_file /= nl) THEN 
   WRITE(ymess,'("Mismatch in # of seaice layers : prep=",I2," nml=",I2)') inl_in_file, nl
   CALL abor1_sfx(ymess)
END IF
 CALL read_surf(&
                     hprogram,'ICENT',int_in_file,iresp)
IF (int_in_file /= nt) THEN
   WRITE(ymess,'("Mismatch in # of seaice categories : prep=",I2," nml=",I2)') int_in_file, nt
   CALL abor1_sfx(ymess)
END IF
!
!*       1.     (Semi-)prognostic fields with only space dimension(s) :
!
 CALL read_surf(&
                     hprogram,'ICEUSTAR',s%TGLT%ust(:,1),iresp)
!
!*       2.     Prognostic fields with space and ice-category dimension(s) :
!
DO jk=1,nt
   WRITE(ylvl,'(I2)') jk
   ycateg='_'//adjustl(ylvl)
   ! .. Read sea ice age for type JK
   CALL read_surf(&
                     hprogram,'ICEAGE'//ycateg,s%TGLT%sit(jk,:,1)%age,iresp)
   ! .. Read melt pond volume for type JK
   CALL read_surf(&
                     hprogram,'ICEVMP'//ycateg,s%TGLT%sit(jk,:,1)%vmp,iresp)
   ! .. Read sea ice surface albedo for type JK
   CALL read_surf(&
                     hprogram,'ICEASN'//ycateg,s%TGLT%sit(jk,:,1)%asn,iresp)
   ! .. Read sea ice fraction for type JK
   CALL read_surf(&
                     hprogram,'ICEFSI'//ycateg, s%TGLT%sit(jk,:,1)%fsi,iresp)
   ! .. Read sea ice thickness for type JK
   CALL read_surf(&
                     hprogram,'ICEHSI'//ycateg, s%TGLT%sit(jk,:,1)%hsi,iresp)
   ! .. Read sea ice salinity for type JK
   CALL read_surf(&
                     hprogram,'ICESSI'//ycateg, s%TGLT%sit(jk,:,1)%ssi,iresp)
   ! .. Read sea ice surface temperature for type JK
   CALL read_surf(&
                     hprogram,'ICETSF'//ycateg, s%TGLT%sit(jk,:,1)%tsf,iresp)
   ! .. Read snow thickness for type JK
   CALL read_surf(&
                     hprogram,'ICEHSN'//ycateg, s%TGLT%sit(jk,:,1)%hsn,iresp)
   ! .. Read snow density for type JK
   CALL read_surf(&
                     hprogram,'ICERSN'//ycateg, s%TGLT%sit(jk,:,1)%rsn,iresp)
   !
   !*       3.     Prognostic fields with space, ice-category and layer dimensions :
   !
   DO jl=1,nl
      WRITE(ylvl,'(I2)') jl
      ylevel=ycateg(1:len_trim(ycateg))//'_'//adjustl(ylvl)   
      ! .. Read sea ice vertical gltools_enthalpy profile for type JK and level JL  
      CALL read_surf(&
                     hprogram,'ICEH'//ylevel, s%TGLT%sil(jl,jk,:,1)%ent,iresp)
   END DO
END DO
!
!    4.  Compute ice class existence boolean from ice fractions:
!
WHERE ( s%TGLT%sit(:,:,1)%fsi<epsil1 )
   s%TGLT%sit(:,:,1)%esi = .false. 
ELSEWHERE 
   s%TGLT%sit(:,:,1)%esi = .true. 
ENDWHERE
!
!    4.1 Run original Gelato checks on values read in restart
!
! .. Detect negative ice concentrations
!
DO jx=1,nx
   DO jl=1,nt 
      IF ( s%TGLT%sit(jl,jx,1)%fsi<0. ) THEN
         WRITE(kluout,*)  &
              '**** WARNING **** Correcting problem in ice conc. < 0 at i=',  &
              1,' j=',jx,' k=',jl
         s%TGLT%sit(jl,jx,1)%fsi = 0.
      ENDIF
   END DO
   !
   zfsit = sum( s%TGLT%sit(:,jx,1)%fsi )
   !
   ! .. Detect total concentrations that exceed unity
   !
   IF ( zfsit>1. ) THEN
      WRITE(kluout,*)  &
           '**** WARNING **** Correcting problem in total ice conc. >1 at i=',  &
           1,' j=',jx,' fsi=',zfsit
      s%TGLT%sit(:,jx,1)%fsi = s%TGLT%sit(:,jx,1)%fsi / zfsit
   ENDIF
   !
   ! .. Detect non zero concentrations but zero thickness (no consequence) 
   !
   WHERE( s%TGLT%sit(:,jx,1)%fsi>epsil1 .AND. s%TGLT%sit(:,jx,1)%hsi<epsil1)
      s%TGLT%sit(:,jx,1)%fsi=0.
      s%TGLT%sit(:,jx,1)%hsi=0.
      s%TGLT%sit(:,jx,1)%hsn=0.
   ENDWHERE
   !
END DO

!    5. Initalize Gelato domain parameters
!
!    All points of Surfex 1D grid in seaflux are sea points
!
s%TGLT%dom(:,1)%tmk=1
s%TGLT%dom(:,1)%imk=1
!
!    Masks for U- and V- grid point are not used
!
s%TGLT%dom(:,1)%umk=1
s%TGLT%dom(:,1)%vmk=1
!
!    lat,lon,srf are inherited from seaflux grid
!
s%TGLT%dom(:,1)%lon=sg%XLON(:)*xpi/180.
s%TGLT%dom(:,1)%lat=sg%XLAT(:)*xpi/180.
!
!    Except in Gelato dynamics, mesh lengths are used only to compute mesh area
!    Hence, a simple setting can be used
!
s%TGLT%dom(:,1)%dxc=sg%XMESH_SIZE(:)**0.5
s%TGLT%dom(:,1)%dyc=s%TGLT%dom(:,1)%dxc
s%TGLT%dom(:,1)%srf=sg%XMESH_SIZE(:)
!
!    Surface of local and global ocean domain (ghost points are masked out)
!
xdomsrf = sum( s%TGLT%dom(:,1)%srf, mask=(s%TGLT%dom(:,1)%tmk==1) )
xdomsrf_g = xdomsrf
#if ! defined in_arpege
 CALL mpp_sum(xdomsrf_g) 
#else
! Avoid zero divide in Gelato computation of global area averages
xdomsrf_g = max(xdomsrf_g, 1.e-9)
#endif
!
!    7. Initalize Gelato time parameters
!
s%TGLT%ind%beg=1
!
!   Dummy high value for end time. Implies only that Gelato won't output 
!   its own format of run-long averaged diagnostics (which are useless 
!   in Surfex diags logic) 
!
s%TGLT%ind%end=50000000 
!
!   8. Initalize Gelato bathymetry - change sign w.r.t Surfex
!
s%TGLT%bat(:,1)=-s%XSEABATHY
!
!
!* Sea ice thickness nudging data
!
IF(s%LINTERPOL_SIT)THEN
   !
   ALLOCATE(s%XFSIT(klu))
   !
   !Precedent, Current, Next, and Second-next Monthly SIT
   inmth=4   
   !
   ALLOCATE(s%XSIT_MTH(klu,inmth))
   DO jmth=1,inmth
      WRITE(ymth,'(I2)') (jmth-1)
      yrecfm='SIT_MTH'//adjustl(ymth(:len_trim(ymth)))
      CALL read_surf(&
                     hprogram,yrecfm,s%XSIT_MTH(:,jmth),iresp)
      CALL check_seaice(yrecfm,s%XSIT_MTH(:,jmth))
   ENDDO
   !
   CALL interpol_sst_mth(s, &
                         s%TTIME%TDATE%YEAR,s%TTIME%TDATE%MONTH,s%TTIME%TDATE%DAY,'H',s%XFSIT)
   !
ELSE
   ! 
   ALLOCATE(s%XFSIT(0))
   ALLOCATE(s%XSIT_MTH(0,0))
   !
ENDIF
!
!! Initialize the coupling variables with 'snapshot' prognostic variables
! (for now, averaged over ice categories)  
!
 CALL glt_sndatmf( s%TGLT, xttsi - xtt )
s%XSIC(:)     = s%TGLT%ice_atm(1,:,1)%fsi 
s%XTICE(:)    = s%TGLT%ice_atm(1,:,1)%tsf 
s%XICE_ALB(:) = s%TGLT%ice_atm(1,:,1)%alb 
!
! Must init ocean mixed layer temp with sensible value for getting correct diag for time step=0
s%TGLT%oce_all(:,1)%tml=s%XSST(:)
!
IF (lhook) CALL dr_hook('READ_SEAICE_n',1,zhook_handle)
!
!-------------------------------------------------------------------------------
 CONTAINS
!-------------------------------------------------------------------------------
!
SUBROUTINE check_seaice(HFIELD,PFIELD)
!
!
IMPLICIT NONE
!
 CHARACTER(LEN=12),  INTENT(IN) :: hfield
REAL, DIMENSION(:), INTENT(IN) :: pfield
!
REAL            :: zmax,zmin
INTEGER         :: ji, ierrc
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('READ_SEAICE_n:CHECK_SEAICE',0,zhook_handle)
!
zmin=-1.0e10
zmax=1.0e10
!
ierrc=0
!
DO ji=1,klu
   IF(pfield(ji)>zmax.OR.pfield(ji)<zmin)THEN
      ierrc=ierrc+1
      WRITE(kluout,*)'PROBLEM FIELD '//trim(hfield)//' =',pfield(ji),&
                     'NOT REALISTIC AT LOCATION (LAT/LON)',sg%XLAT(ji),sg%XLON(ji)
   ENDIF
ENDDO
!         
IF(ierrc>0) CALL abor1_sfx('READ_SEAICE_n: FIELD '//trim(hfield)//' NOT REALISTIC')
!
IF (lhook) CALL dr_hook('READ_SEAICE_n:CHECK_SEAICE',1,zhook_handle)

END SUBROUTINE check_seaice
!
!------------------------------------------------------------------------------
!
END SUBROUTINE read_seaice_n