surface_fields.F90

Go to the documentation of this file.

MODULE surface_fields

!     Purpose.
!     --------

!      SURFACE_FIELDS contains data structures and manipulation routines
!      for the surface (physics) fields in the IFS            

!      This module is a mix of declarations, type definitions and 
!       subroutines linked with surface fields. There are four parts:
!       1/ Declaration of dimensions (including some parameter variables).
!       2/ Definition of types.
!       3/ Declarations:
!          Declaration of variables SP_[group], YSP_[group]D, YSP_[group]
!           (prognostic surface fields).
!          Declaration of variables SD_[group], YSD_[group]D, YSD_[group]
!           (diagnostic surface fields).
!       4/ Some routines linked to the surface data flow:
!          * INI_SFLP3: Initialize 3-D surface field group
!          * SETUP_SFLP3: Setup 3-D surface field
!          * INI_SFLP2: Initialize 2-D surface field group
!          * SETUP_SFLP2: Setup 2-D surface field
!          * GPPOPER: Operations on prognostic surface fields
!          * GPOPER: Operations on ALL surface groups
!          * GPOPER_2: Operations on 2-D surface groups
!          * GPOPER_3: Operations on 3-D surface groups
!          * SURF_STORE: Store all surface fields
!          * SURF_RESTORE: Restore all surface fields
!          * ALLO_SURF: Allocate surface field arrays
!          * DEALLO_SURF: Deallocate surface field arrays

!     Author.
!     -------
!     Mats Hamrud(ECMWF)

!     Modifications.
!     --------------
!        Original : 2006-07-01
!        Modifications:
!        K. Yessad (25 Oct 2006): rephase ALARO0 contribution.
!        K. Yessad (26 Oct 2006): add missing comments.

!-------------------------------------------------------------------------

USE parkind1  ,ONLY : jpim     ,jprb
USE yomdim    ,ONLY : nproma, ngpblks, nundefld
USE yomlun    ,ONLY : nulout, nulerr
USE yomct0    ,ONLY : ltwotl
USE yomdyn   , ONLY : repsp1
USE yomhook   ,ONLY : lhook,   dr_hook
IMPLICIT NONE
SAVE

!#include "abor1.intfb.h"
!     -------------------------------------------------------------------------

INTEGER(KIND=JPIM), PARAMETER :: jpmaxsflds=100 ! Max number of fields in individual group
INTEGER(KIND=JPIM), PARAMETER :: jpmaxstraj=100 ! Dimension of NSTRAJGRIB
INTEGER(KIND=JPIM) :: nsurf=0              ! Number of surf var.
INTEGER(KIND=JPIM) :: nsurfl=0             ! Number of surf flds (fields*levels)
INTEGER(KIND=JPIM) :: ndimsurf=0           ! Total of surf var (includes timelevels etc)
INTEGER(KIND=JPIM) :: ndimsurfl=0          ! Total dimension of all surface variables
INTEGER(KIND=JPIM) :: nprogsurf=0          ! Number of prognostic surf var.
INTEGER(KIND=JPIM) :: nprogsurfl=0         ! Number of prognostic surf flds (fields*levels)
INTEGER(KIND=JPIM) :: nofftraj             ! Offset in surf trajectory
INTEGER(KIND=JPIM) :: nofftraj_cst         ! Offset in "constant" surf trajectory
INTEGER(KIND=JPIM) :: nptrsurf             ! Used by routine GPOPER
INTEGER(KIND=JPIM) :: nstrajgrib(jpmaxstraj) ! Used in trajectory setup

REAL(KIND=JPRB),ALLOCATABLE   :: surf_store_array(:,:,:) ! Backup array for surf (see routineSURF_STORE )
! General type defintions

! 2D surface field structure
TYPE type_surf_mtl_2d
INTEGER(KIND=JPIM) :: mp                   ! Basic field pointer
INTEGER(KIND=JPIM) :: mp0                  ! Field pointer timelevel  0 (prognostic fields)
INTEGER(KIND=JPIM) :: mp9                  ! Field pointer timelevel -1 (prognostic fields)
INTEGER(KIND=JPIM) :: mp1                  ! Field pointer timelevel +1 (prognostic fields)
INTEGER(KIND=JPIM) :: mp5                  ! Field pointer trajectory
INTEGER(KIND=JPIM) :: igrbcode             ! GRIB parameter code (default: -999)
CHARACTER(LEN=16)  :: cname                ! ARPEGE field name   (default: all spaces)
REAL(KIND=JPRB)    :: refvali              ! Default value       (default: 0.0)
INTEGER(KIND=JPIM) :: nreqin               ! -1 - initial value from default (default)
                                           ! +1 - initial value from reading file
                                           !  0 - no initial value
INTEGER(KIND=JPIM) :: itraj                !  0 not in trajectory (default)
                                           !  1 in trajectory
                                           !  2 in "constant" trajectory
END TYPE type_surf_mtl_2d

! 3D surface field structure
TYPE type_surf_mtl_3d
INTEGER(KIND=JPIM) :: mp   ! Basic field pointer
INTEGER(KIND=JPIM) :: mp0  ! Field pointer timelevel  0 (prognostic fields)
INTEGER(KIND=JPIM) :: mp9  ! Field pointer timelevel -1 (prognostic fields)
INTEGER(KIND=JPIM) :: mp1  ! Field pointer timelevel +1 (prognostic fields)
INTEGER(KIND=JPIM) :: mp5  ! Field pointer trajectory
INTEGER(KIND=JPIM),POINTER :: igrbcode(:)  ! GRIB parameter code (default: -999)
CHARACTER(LEN=16) ,POINTER :: cname(:)     ! ARPEGE field name   (default: all spaces)
REAL(KIND=JPRB)   ,POINTER :: refvali(:)   ! Default value       (default: 0.0)
INTEGER(KIND=JPIM),POINTER :: nreqin(:)    ! -1 - initial value from default (default)
                                           ! +1 - initial value from reading file
                                           !  0 - no initial value
INTEGER(KIND=JPIM) :: itraj                !  0 not in trajectory (default)
                                           !  1 in trajectory
                                           !  2 in "constant" trajectory
END TYPE type_surf_mtl_3d

! Descriptor pertaining to group
TYPE type_surf_gen
INTEGER(KIND=JPIM) :: numflds         ! Number of field in group
INTEGER(KIND=JPIM) :: ndim            ! Field dimenion
INTEGER(KIND=JPIM) :: nlevs           ! Number of levels (for multi level groups)
INTEGER(KIND=JPIM) :: iptr            ! Internal use
INTEGER(KIND=JPIM) :: iptr5           ! Internal use
INTEGER(KIND=JPIM) :: ndim5           ! Dimension of trajectory array
INTEGER(KIND=JPIM) :: nofftraj        ! Internal use
INTEGER(KIND=JPIM) :: nofftraj_cst    ! Internal use
CHARACTER(LEN=16)  :: cgrpname        ! Name of group (for prints)
LOGICAL            :: l3d             ! TRUE if multi-level field (3-D)
LOGICAL            :: lmtl            ! TRUE if prognostic field (multi time level)
END TYPE type_surf_gen

! Type descriptor for derived type for communicating with GPOPER (see below)
TYPE type_sfl_comm
INTEGER(KIND=JPIM) :: igrbcode
LOGICAL            :: l_ok 
CHARACTER(LEN=16)  :: cname
INTEGER(KIND=JPIM) :: ifldnum
REAL(KIND=JPRB)    :: value
INTEGER(KIND=JPIM) :: iptrsurf
INTEGER(KIND=JPIM) :: icodes(jpmaxsflds)
INTEGER(KIND=JPIM) :: icount
END TYPE type_sfl_comm

! Group specific type definitions

! * Group SB=SOILB: soil prognostic quantities for the different reservoirs
!    (four reservoirs at ECMWF, deep reservoir at METEO-FRANCE):
TYPE type_sfl_soilb
TYPE(type_surf_mtl_3d),POINTER :: yt    ! temperature
TYPE(type_surf_mtl_3d),POINTER :: yq    ! liquid water content
TYPE(type_surf_mtl_3d),POINTER :: ytl   ! ice water content (for MF)
TYPE(type_surf_mtl_3d),POINTER :: ysb(:)
END TYPE type_sfl_soilb

! * Group SG=SNOWG: surface snow prognostic quantities:
TYPE type_sfl_snowg
TYPE(type_surf_mtl_2d),POINTER :: yf    ! content of surface snow
TYPE(type_surf_mtl_2d),POINTER :: ya    ! snow albedo
TYPE(type_surf_mtl_2d),POINTER :: yr    ! snow density
TYPE(type_surf_mtl_2d),POINTER :: yt    ! total albedo (diagnostic for MF for LVGSN)
TYPE(type_surf_mtl_2d),POINTER :: ysg(:)
END TYPE type_sfl_snowg

! * Group RR=RESVR: surface prognostic quantities (ECMWF) or
!   surface + superficial reservoir prognostic quantities (MF):
!   Remark:
!    at ECMWF there are 4 soil reservoirs and there is a
!    clear distinction between the soil reservoirs (group SOILB)
!    and the surface (group RESVR);
!    at METEO-FRANCE there is a deep reservoir (group SOILB) and a
!    superficial reservoir (group RESVR):
!    - there is a skin surface temperature (Ts) which is the temperature at the
!      interface surface/superficial reservoir (and not two separate quantities
!      for superficial reservoir and surface)
!    - there is a skin surface water content (denoted by Wl) and a superficial
!      reservoir water content (denoted by Ws).
!    - there is a superficial reservoir ice content but no surface ice content.
!    (remark k.y.: it would have been more logical to use group name
!    RESVR for internal reservoirs and group name SOILB for surface!).
TYPE type_sfl_resvr
TYPE(type_surf_mtl_2d),POINTER :: yt    ! skin temperature (Ts)
TYPE(type_surf_mtl_2d),POINTER :: yw    ! skin water content (Wskin) at ECMWF
                                        ! superficial reservoir water content (Ws) at MF
TYPE(type_surf_mtl_2d),POINTER :: yfc   ! skin water content (Wl) at MF
TYPE(type_surf_mtl_2d),POINTER :: yic   ! superficial reservoir ice
TYPE(type_surf_mtl_2d),POINTER :: yfp1  ! interpolated Ts for 2nd part of 927-FULLPOS
TYPE(type_surf_mtl_2d),POINTER :: yrr(:)
END TYPE type_sfl_resvr

! * Group WS=WAVES: surface prognostic quantities over sea:
TYPE type_sfl_waves
TYPE(type_surf_mtl_2d),POINTER :: yws(:)
TYPE(type_surf_mtl_2d),POINTER :: ychar ! Charnock constant
END TYPE type_sfl_waves

! * Group EP=EXTRP: extra 3-d prognostic fields:
TYPE type_sfl_extrp
TYPE(type_surf_mtl_3d),POINTER :: yep(:)
END TYPE type_sfl_extrp

! * Group X2=XTRP2: extra 2-d prognostic fields:
!   (is used for precipitation fields in CANARI)
TYPE type_sfl_xtrp2
TYPE(type_surf_mtl_2d),POINTER :: yx2(:)
END TYPE type_sfl_xtrp2

! * Group CI=CANRI: 2-d prognostic fields for CANARI:
TYPE type_sfl_canri
TYPE(type_surf_mtl_2d),POINTER :: yci(:)
END TYPE type_sfl_canri

! * Group VF=VARSF: climatological/geographical diagnostic fields:
TYPE type_sfl_varsf
TYPE(type_surf_mtl_2d),POINTER :: yz0f    ! gravity * surface roughness length
TYPE(type_surf_mtl_2d),POINTER :: yalbf   ! surface shortwave albedo
TYPE(type_surf_mtl_2d),POINTER :: yemisf  ! surface longwave emissivity
TYPE(type_surf_mtl_2d),POINTER :: ygetrl  ! standard deviation of orography
TYPE(type_surf_mtl_2d),POINTER :: yitm    ! land-sea mask
TYPE(type_surf_mtl_2d),POINTER :: yveg    ! vegetation cover
TYPE(type_surf_mtl_2d),POINTER :: yvrlan  ! anisotropy of the sub-grid scale orography
TYPE(type_surf_mtl_2d),POINTER :: yvrldi  ! angle of the direction of orography with the x axis
TYPE(type_surf_mtl_2d),POINTER :: ysig    ! characteristic orographic slope
TYPE(type_surf_mtl_2d),POINTER :: yalbsf  ! soil shortwave albedo
TYPE(type_surf_mtl_2d),POINTER :: ycont   ! fraction of land
TYPE(type_surf_mtl_2d),POINTER :: ysst    ! (open) sea surface temperature
TYPE(type_surf_mtl_2d),POINTER :: ylz0h   ! logarithm of roughness length for heat
TYPE(type_surf_mtl_2d),POINTER :: ycvl    ! low vegetation cover
TYPE(type_surf_mtl_2d),POINTER :: ycvh    ! high vegetation cover
TYPE(type_surf_mtl_2d),POINTER :: ytvl    ! low vegetation type
TYPE(type_surf_mtl_2d),POINTER :: ytvh    ! high vegetation type
TYPE(type_surf_mtl_2d),POINTER :: yci     ! sea ice fraction
TYPE(type_surf_mtl_2d),POINTER :: yucur   ! U-component of the ocean current
TYPE(type_surf_mtl_2d),POINTER :: yvcur   ! V-component of the ocean current
TYPE(type_surf_mtl_2d),POINTER :: yz0rlf  ! gravity * vegetation roughness length
TYPE(type_surf_mtl_2d),POINTER :: yco2o   ! oceanic CO2 flux
TYPE(type_surf_mtl_2d),POINTER :: yco2b   ! biosphere CO2 flux
TYPE(type_surf_mtl_2d),POINTER :: yco2a   ! anthropogenic CO2 flux
TYPE(type_surf_mtl_2d),POINTER :: ysdfor  ! SD filtered orography
TYPE(type_surf_mtl_2d),POINTER :: yaluvp  ! MODIS-derived parallel albedo for shortwave radiation
TYPE(type_surf_mtl_2d),POINTER :: yaluvd  ! MODIS-derived diffuse albedo for shortwave radiation
TYPE(type_surf_mtl_2d),POINTER :: yalnip  ! MODIS-derived parallel albedo for longwave radiation
TYPE(type_surf_mtl_2d),POINTER :: yalnid  ! MODIS-derived diffuse albedo for longwave radiation
TYPE(type_surf_mtl_2d),POINTER :: ysf6    ! anthropogenic SF6 flux
TYPE(type_surf_mtl_2d),POINTER :: yfp1    ! surface orography in the 2nd part of FULLPOS-927
TYPE(type_surf_mtl_2d),POINTER :: yvf(:)
END TYPE type_sfl_varsf

! * Group VP=VCLIP: deep soil diagnostic fields:
TYPE type_sfl_vclip
TYPE(type_surf_mtl_2d),POINTER :: ytpc    ! climatological deep layer temperature
TYPE(type_surf_mtl_2d),POINTER :: ywpc    ! climatological deep layer moisture
TYPE(type_surf_mtl_2d),POINTER :: yvp(:)
END TYPE type_sfl_vclip

! * Group VV=VCLIV: vegetation diagnostic fields:
TYPE type_sfl_vcliv
TYPE(type_surf_mtl_2d),POINTER :: yarg    ! silt percentage within soil
TYPE(type_surf_mtl_2d),POINTER :: ysab    ! percentage of sand within the soil
TYPE(type_surf_mtl_2d),POINTER :: yd2     ! soil depth
TYPE(type_surf_mtl_2d),POINTER :: yiveg   ! type of vegetation
TYPE(type_surf_mtl_2d),POINTER :: yrsmin  ! stomatal minimum resistance
TYPE(type_surf_mtl_2d),POINTER :: ylai    ! leaf area index
TYPE(type_surf_mtl_2d),POINTER :: yhv     ! resistance to evapotranspiration
TYPE(type_surf_mtl_2d),POINTER :: yz0h    ! gravity * roughness length for heat
TYPE(type_surf_mtl_2d),POINTER :: yals    ! albedo of bare ground
TYPE(type_surf_mtl_2d),POINTER :: yalv    ! albedo of vegetation
TYPE(type_surf_mtl_2d),POINTER :: yvv(:)
END TYPE type_sfl_vcliv

! * Group VN=VCLIN: cloudiness diagnostic predictors:
TYPE type_sfl_vclin
TYPE(type_surf_mtl_2d),POINTER :: ytop    ! index of convective cloud top
TYPE(type_surf_mtl_2d),POINTER :: ybas    ! index of convective cloud base
TYPE(type_surf_mtl_2d),POINTER :: yacpr   ! averaged convective precipitaion rate
TYPE(type_surf_mtl_2d),POINTER :: yvn(:)
END TYPE type_sfl_vclin

! * Group VH=VCLIH: convective cloud diagnostic fields:
TYPE type_sfl_vclih
TYPE(type_surf_mtl_2d),POINTER :: ytcch   ! total convective cloudiness
TYPE(type_surf_mtl_2d),POINTER :: yscch   ! convective cloud summit
TYPE(type_surf_mtl_2d),POINTER :: ybcch   ! convective cloud base
TYPE(type_surf_mtl_2d),POINTER :: ypblh   ! PBL height
TYPE(type_surf_mtl_2d),POINTER :: yspsh   ! variable for prognostic convection scheme (ALARO)
TYPE(type_surf_mtl_2d),POINTER :: yvh(:)
END TYPE type_sfl_vclih

! * Group VA=VCLIA: aerosol diagnostic fields:
TYPE type_sfl_vclia
TYPE(type_surf_mtl_2d),POINTER :: ysea    ! aerosol: sea
TYPE(type_surf_mtl_2d),POINTER :: ylan    ! aerosol: land
TYPE(type_surf_mtl_2d),POINTER :: ysoo    ! aerosol: soot
TYPE(type_surf_mtl_2d),POINTER :: ydes    ! aerosol: desert
TYPE(type_surf_mtl_2d),POINTER :: ysul    ! aerosol: sulfate
TYPE(type_surf_mtl_2d),POINTER :: yvol    ! aerosol: volcano
TYPE(type_surf_mtl_2d),POINTER :: ynud    ! aerosol: nudging
TYPE(type_surf_mtl_2d),POINTER :: yva(:)
END TYPE type_sfl_vclia

! * Group VG=VCLIG: ice-coupler diagnostic fields:
TYPE type_sfl_vclig
TYPE(type_surf_mtl_2d),POINTER :: yicfr   ! sea-ice fraction
TYPE(type_surf_mtl_2d),POINTER :: ysoup   ! upward solar flux over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: yirup   ! upward IR flux over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: ychss   ! sensible heat over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: yevap   ! evaporation over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: ytaux   ! U-component of stress over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: ytauy   ! V-component of stress over sea-ice
TYPE(type_surf_mtl_2d),POINTER :: yvg(:)
END TYPE type_sfl_vclig

! * Group VC=VO3ABC: A,B and C (Climatological ozone profiles) diagnostic fields:
TYPE type_sfl_vo3abc
TYPE(type_surf_mtl_2d),POINTER :: ya      ! A climatological ozone profile
TYPE(type_surf_mtl_2d),POINTER :: yb      ! B climatological ozone profile
TYPE(type_surf_mtl_2d),POINTER :: yc      ! C climatological ozone profile
TYPE(type_surf_mtl_2d),POINTER :: yvc(:)
END TYPE type_sfl_vo3abc

! * Group VD=VDIAG: (ECMWF) diagnostic fields:
TYPE type_sfl_vdiag
TYPE(type_surf_mtl_2d),POINTER :: ylsp   !Large scale precipitation
TYPE(type_surf_mtl_2d),POINTER :: ycp    !Convective precipitation
TYPE(type_surf_mtl_2d),POINTER :: ysf    !Snowfall
TYPE(type_surf_mtl_2d),POINTER :: ybld   !Boundary layer dissipation
TYPE(type_surf_mtl_2d),POINTER :: ysshf  !Surface sensible heat flux
TYPE(type_surf_mtl_2d),POINTER :: yslhf  !Surface latent heat flux
TYPE(type_surf_mtl_2d),POINTER :: ymsl   !Mean sea level pressure
TYPE(type_surf_mtl_2d),POINTER :: ytcc   !Total cloud cover
TYPE(type_surf_mtl_2d),POINTER :: y10u   !U-wind at 10 m
TYPE(type_surf_mtl_2d),POINTER :: y10v   !V-wind at 10 m
TYPE(type_surf_mtl_2d),POINTER :: y2t    !Temperature at 2 m
TYPE(type_surf_mtl_2d),POINTER :: y2d    !Dewpoint temperature at 2 m
TYPE(type_surf_mtl_2d),POINTER :: yssr   !Surface solar radiation
TYPE(type_surf_mtl_2d),POINTER :: ystr   !Surface thermal radiation
TYPE(type_surf_mtl_2d),POINTER :: ytsr   !Top solar radiation
TYPE(type_surf_mtl_2d),POINTER :: yttr   !Top thermal radiation
TYPE(type_surf_mtl_2d),POINTER :: yewss  !Instantaneous surface U-wind stress
TYPE(type_surf_mtl_2d),POINTER :: ynsss  !Instantaneous surface V-wind stress
TYPE(type_surf_mtl_2d),POINTER :: ye     !Water evaporation
TYPE(type_surf_mtl_2d),POINTER :: yccc   !Convective cloud cover
TYPE(type_surf_mtl_2d),POINTER :: ylcc   !Low cloud cover
TYPE(type_surf_mtl_2d),POINTER :: ymcc   !Medium cloud cover
TYPE(type_surf_mtl_2d),POINTER :: yhcc   !High cloud cover
TYPE(type_surf_mtl_2d),POINTER :: ylgws  !Zonal gravity wave stress
TYPE(type_surf_mtl_2d),POINTER :: ymgws  !Meridian gravity wave stress
TYPE(type_surf_mtl_2d),POINTER :: ygwd   !Gravity wave dissipation
TYPE(type_surf_mtl_2d),POINTER :: ymx2t  !Maximum temperature at 2 m
TYPE(type_surf_mtl_2d),POINTER :: ymn2t  !Minimum temperature at 2 m
TYPE(type_surf_mtl_2d),POINTER :: yro    !Runoff
TYPE(type_surf_mtl_2d),POINTER :: yalb   !(surface shortwave) albedo
TYPE(type_surf_mtl_2d),POINTER :: yiewss !Instantaneous surface zonal component of stress
TYPE(type_surf_mtl_2d),POINTER :: yinsss !Instantaneous surface meridian component of stress
TYPE(type_surf_mtl_2d),POINTER :: yisshf !Instantaneous surface heat flux
TYPE(type_surf_mtl_2d),POINTER :: yie    !Instantaneous surface moisture flux
TYPE(type_surf_mtl_2d),POINTER :: ycsf   !Convective snow fall
TYPE(type_surf_mtl_2d),POINTER :: ylssf  !Large scale snowfall
TYPE(type_surf_mtl_2d),POINTER :: yz0f   !Gravity * surface roughness length
TYPE(type_surf_mtl_2d),POINTER :: ylz0h  !Logarithm of z0 times heat flux
TYPE(type_surf_mtl_2d),POINTER :: ytcw   !Total water content in a vertical column
TYPE(type_surf_mtl_2d),POINTER :: ytcwv  !Total water vapor content in a vertical column
TYPE(type_surf_mtl_2d),POINTER :: ytclw  !Total liquid water content in a vertical column
TYPE(type_surf_mtl_2d),POINTER :: ytciw  !Total ice water content in a vertical column
TYPE(type_surf_mtl_2d),POINTER :: yssrd  !Downward surface solar radiation
TYPE(type_surf_mtl_2d),POINTER :: ystrd  !Downward surface thermic radiation
TYPE(type_surf_mtl_2d),POINTER :: yblh   !Height of boundary layer
TYPE(type_surf_mtl_2d),POINTER :: ysund  !Sunshine duration
TYPE(type_surf_mtl_2d),POINTER :: yspar  !Surface downward PARadiation
TYPE(type_surf_mtl_2d),POINTER :: ysuvb  !Surface downward UV-B radiation
TYPE(type_surf_mtl_2d),POINTER :: ycape  !Conv.avail.potential energy (CAPE)
TYPE(type_surf_mtl_2d),POINTER :: ytsrc  !Top solar radiation clear sky
TYPE(type_surf_mtl_2d),POINTER :: yttrc  !Top thermal radiation clear sky
TYPE(type_surf_mtl_2d),POINTER :: yssrc  !Surface solar radiation clear sky
TYPE(type_surf_mtl_2d),POINTER :: ystrc  !Surface thermal radiation clear sky
TYPE(type_surf_mtl_2d),POINTER :: yes    !Evaporation of snow
TYPE(type_surf_mtl_2d),POINTER :: ysmlt  !Snow melt
TYPE(type_surf_mtl_2d),POINTER :: y10fg  !Wind gust at 10 m (max since previous pp)
TYPE(type_surf_mtl_2d),POINTER :: ylspf  !Large scale precipitation fraction
TYPE(type_surf_mtl_2d),POINTER :: ytco3  !Total ozone content in a vertical column
TYPE(type_surf_mtl_2d),POINTER :: yvimd  !Vertically integrated mass divergence
TYPE(type_surf_mtl_2d),POINTER :: ysparc !Surface clear-sky parallel radiation
TYPE(type_surf_mtl_2d),POINTER :: ystinc !TOA (top of atmosph?) incident solar radiation
TYPE(type_surf_mtl_2d),POINTER :: ytcghg(:)  !Total column greenhouse gases
TYPE(type_surf_mtl_2d),POINTER :: ytcgrg(:)  !Total column reactive gases
TYPE(type_surf_mtl_2d),POINTER :: ytctrac(:) !Total column tracers
TYPE(type_surf_mtl_2d),POINTER :: yvd(:)
END TYPE type_sfl_vdiag

! * Group VX=VCLIX: auxilary climatological diagnostic fields:
TYPE type_sfl_vclix
TYPE(type_surf_mtl_2d),POINTER :: yoro    ! climatological surface geopotential
TYPE(type_surf_mtl_2d),POINTER :: ytsc    ! climatological surface temperature
TYPE(type_surf_mtl_2d),POINTER :: ypws    ! climatological surface max. prop. moisture
TYPE(type_surf_mtl_2d),POINTER :: ypwp    ! climatological deep soil max. prop. moisture
TYPE(type_surf_mtl_2d),POINTER :: ysno    ! climatological snow cover
TYPE(type_surf_mtl_2d),POINTER :: ytpc    ! climatological deep soil temperature
TYPE(type_surf_mtl_2d),POINTER :: ysab    ! climatologic percentage of sand within the soil
TYPE(type_surf_mtl_2d),POINTER :: yxd2    ! climatologic soil depth
TYPE(type_surf_mtl_2d),POINTER :: ylsm    ! climatologic land sea mask
TYPE(type_surf_mtl_2d),POINTER :: yiveg   ! climatologic type of vegetation
TYPE(type_surf_mtl_2d),POINTER :: yvx(:)
END TYPE type_sfl_vclix

! * Group XA=VEXTRA: extra 3-d diagnostic fields:
TYPE type_sfl_vextra
TYPE(type_surf_mtl_3d),POINTER :: yxa(:)
END TYPE type_sfl_vextra

! * Group X2=VEXTR2: extra 2-d diagnostic fields:
TYPE type_sfl_vextr2
TYPE(type_surf_mtl_2d),POINTER :: yx2(:)
END TYPE type_sfl_vextr2

! End of type definitions

! Data structures

! Prognostic (multi time level) fields

! Soilb
REAL(KIND=JPRB),ALLOCATABLE :: sp_sb (:,:,:,:)
TYPE(type_surf_gen)    :: ysp_sbd
TYPE(type_sfl_soilb)   :: ysp_sb

! Snowg
REAL(KIND=JPRB),ALLOCATABLE :: sp_sg (:,:,:)
TYPE(type_surf_gen)    :: ysp_sgd
TYPE(type_sfl_snowg)   :: ysp_sg

! Resvr
REAL(KIND=JPRB),ALLOCATABLE :: sp_rr (:,:,:)
TYPE(type_surf_gen)    :: ysp_rrd
TYPE(type_sfl_resvr)   :: ysp_rr


! Extrp
REAL(KIND=JPRB),ALLOCATABLE :: sp_ep (:,:,:,:)
TYPE(type_surf_gen)    :: ysp_epd
TYPE(type_sfl_extrp)   :: ysp_ep 

! Xtrp2
REAL(KIND=JPRB),ALLOCATABLE :: sp_x2 (:,:,:)
TYPE(type_surf_gen)    :: ysp_x2d
TYPE(type_sfl_xtrp2)   :: ysp_x2

! Canri
REAL(KIND=JPRB),ALLOCATABLE :: sp_ci (:,:,:)
TYPE(type_surf_gen)    :: ysp_cid
TYPE(type_sfl_canri)   :: ysp_ci

! One time level fields

! Varsf
REAL(KIND=JPRB),ALLOCATABLE :: sd_vf (:,:,:)
TYPE(type_surf_gen)    :: ysd_vfd
TYPE(type_sfl_varsf)   :: ysd_vf

! Vclip
REAL(KIND=JPRB),ALLOCATABLE :: sd_vp (:,:,:)
TYPE(type_surf_gen)    :: ysd_vpd
TYPE(type_sfl_vclip)   :: ysd_vp

! Vcliv
REAL(KIND=JPRB),ALLOCATABLE :: sd_vv (:,:,:)
TYPE(type_surf_gen)    :: ysd_vvd
TYPE(type_sfl_vcliv)   :: ysd_vv

! Vclin
REAL(KIND=JPRB),ALLOCATABLE :: sd_vn (:,:,:)
TYPE(type_surf_gen)    :: ysd_vnd
TYPE(type_sfl_vclin)   :: ysd_vn

! Vclih
REAL(KIND=JPRB),ALLOCATABLE :: sd_vh (:,:,:)
TYPE(type_surf_gen)    :: ysd_vhd
TYPE(type_sfl_vclih)   :: ysd_vh

! Vclia
REAL(KIND=JPRB),ALLOCATABLE :: sd_va (:,:,:)
TYPE(type_surf_gen)    :: ysd_vad
TYPE(type_sfl_vclia)   :: ysd_va

! Vo3abc
REAL(KIND=JPRB),ALLOCATABLE :: sd_vc (:,:,:)
TYPE(type_surf_gen)    :: ysd_vcd
TYPE(type_sfl_vo3abc)  :: ysd_vc

! Vdiag
REAL(KIND=JPRB),ALLOCATABLE :: sd_vd (:,:,:)
TYPE(type_surf_gen)    :: ysd_vdd
TYPE(type_sfl_vdiag)   :: ysd_vd

! Waves
REAL(KIND=JPRB),ALLOCATABLE :: sd_ws (:,:,:)
TYPE(type_surf_gen)    :: ysd_wsd
TYPE(type_sfl_waves)   :: ysd_ws

! Vclix
REAL(KIND=JPRB),ALLOCATABLE :: sd_vx (:,:,:)
TYPE(type_surf_gen)    :: ysd_vxd
TYPE(type_sfl_vclix)   :: ysd_vx

! Vextra

REAL(KIND=JPRB),ALLOCATABLE :: sd_xa (:,:,:,:)
TYPE(type_surf_gen)    :: ysd_xad
TYPE(type_sfl_vextra)  :: ysd_xa

! Vextr2

REAL(KIND=JPRB),ALLOCATABLE :: sd_x2 (:,:,:)
TYPE(type_surf_gen)    :: ysd_x2d
TYPE(type_sfl_vextr2)  :: ysd_x2

!$OMP THREADPRIVATE(ndimsurf,ndimsurfl,nofftraj,nofftraj_cst,nprogsurf)
!$OMP THREADPRIVATE(nprogsurfl,nptrsurf,nstrajgrib,nsurf,nsurfl,ysd_va,ysd_vad)
!$OMP THREADPRIVATE(ysd_vc,ysd_vcd,ysd_vd,ysd_vdd,ysd_vf,ysd_vfd,ysd_vh,ysd_vhd)
!$OMP THREADPRIVATE(ysd_vn,ysd_vnd,ysd_vp,ysd_vpd,ysd_vv,ysd_vvd,ysd_vx,ysd_vxd)
!$OMP THREADPRIVATE(ysd_ws,ysd_wsd,ysd_x2,ysd_x2d,ysd_xa,ysd_xad,ysp_ci,ysp_cid)
!$OMP THREADPRIVATE(ysp_ep,ysp_epd,ysp_rr,ysp_rrd,ysp_sb,ysp_sbd,ysp_sg,ysp_sgd)
!$OMP THREADPRIVATE(ysp_x2,ysp_x2d)

!$OMP THREADPRIVATE(sd_va,sd_vc,sd_vd,sd_vf,sd_vh,sd_vn,sd_vp,sd_vv,sd_vx,sd_ws)
!$OMP THREADPRIVATE(sd_x2,sd_xa,sp_ci,sp_ep,sp_rr,sp_sb,sp_sg,sp_x2,surf_store_array)


!-------------------------------------------------------------------------

CONTAINS

!=========================================================================

SUBROUTINE ini_sflp3(YDSC,YD,KFLDS,KLEVS,LDMTL,CDGRPNAME)
! Initialize 3-D surface field group
TYPE(type_surf_gen),INTENT(INOUT)    :: YDSC
TYPE(type_surf_mtl_3d),INTENT(INOUT) :: YD(:)
INTEGER(KIND=JPIM),INTENT(IN)        :: KFLDS
INTEGER(KIND=JPIM),INTENT(IN)        :: KLEVS
LOGICAL,INTENT(IN)                   :: LDMTL
CHARACTER(LEN=*),INTENT(IN)          :: CDGRPNAME

INTEGER(KIND=JPIM) :: JFLD, IMAXF
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:INI_SFLP3',0,zhook_handle)

imaxf = SIZE(yd)
ydsc%NUMFLDS = kflds
ydsc%NLEVS = klevs
ydsc%IPTR  = 1
ydsc%LMTL  = ldmtl
ydsc%CGRPNAME = cdgrpname
ydsc%NDIM5 = 0
ydsc%NOFFTRAJ = nofftraj
ydsc%NOFFTRAJ_CST = nofftraj_cst

nsurf = nsurf+ydsc%NUMFLDS
nsurfl = nsurfl+ydsc%NUMFLDS*ydsc%NLEVS
IF(ldmtl) THEN
  nprogsurf = nprogsurf+ydsc%NUMFLDS
  nprogsurfl = nprogsurfl+ydsc%NUMFLDS*ydsc%NLEVS
ENDIF

IF(ldmtl) THEN
  IF (ltwotl) THEN
    ydsc%NDIM = 2*ydsc%NUMFLDS
  ELSE
    ydsc%NDIM = 3*ydsc%NUMFLDS
  ENDIF
ELSE
  ydsc%NDIM = ydsc%NUMFLDS
ENDIF  
ndimsurf = ndimsurf + ydsc%NDIM    
ndimsurfl = ndimsurfl + ydsc%NDIM*ydsc%NLEVS
    
DO jfld=1,kflds
  ALLOCATE(yd(jfld)%IGRBCODE(klevs))
  ALLOCATE(yd(jfld)%CNAME(klevs))
  ALLOCATE(yd(jfld)%REFVALI(klevs))
  ALLOCATE(yd(jfld)%NREQIN(klevs))
  yd(jfld)%IGRBCODE(:) = -999
  yd(jfld)%CNAME(:) = ''
  yd(jfld)%REFVALI(:) = 0.0_jprb
  yd(jfld)%NREQIN(:) = -1
  yd(jfld)%MP = jfld
  IF (ydsc%LMTL) THEN
    yd(jfld)%MP0 = yd(jfld)%MP
    IF(ltwotl) THEN
      yd(jfld)%MP9 = yd(jfld)%MP0
      yd(jfld)%MP1 = yd(jfld)%MP0+ydsc%NUMFLDS
    ELSE
      yd(jfld)%MP9 = yd(jfld)%MP0+ydsc%NUMFLDS
      yd(jfld)%MP1 = yd(jfld)%MP0+2*ydsc%NUMFLDS
    ENDIF
  ELSE
    yd(jfld)%MP0 = nundefld
    yd(jfld)%MP9 = nundefld
    yd(jfld)%MP1 = nundefld
  ENDIF
  yd(jfld)%MP5 = nundefld
  yd(jfld)%ITRAJ = 0
ENDDO

DO jfld=kflds+1,imaxf
  yd(jfld)%MP  = nundefld
  yd(jfld)%MP0 = nundefld
  yd(jfld)%MP9 = nundefld
  yd(jfld)%MP1 = nundefld
  yd(jfld)%MP5 = nundefld
  yd(jfld)%ITRAJ = 0
ENDDO

WRITE(nulout,*) 'INITIALIZING 3-D SURFACE FIELD GROUP ', ydsc%CGRPNAME
WRITE(nulout,*) 'NUMFLDS=',ydsc%NUMFLDS,' NLEVS=',ydsc%NLEVS,' LMTL=',ydsc%LMTL 

IF (lhook) CALL dr_hook('SURFACE_FIELDS:INI_SFLP3',1,zhook_handle)
END SUBROUTINE ini_sflp3

!=========================================================================

SUBROUTINE setup_sflp3(YDSC,YD,KGRIB,CDNAME,PDEFAULT,KTRAJ,KREQIN)
! Setup 3-D surface field 
TYPE(type_surf_gen),INTENT(INOUT)      :: YDSC
TYPE(type_surf_mtl_3d),INTENT(INOUT)   :: YD
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KGRIB(:)
CHARACTER(LEN=16) ,OPTIONAL,INTENT(IN) :: CDNAME(:)
REAL(KIND=JPRB)   ,OPTIONAL,INTENT(IN) :: PDEFAULT(:)
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KTRAJ
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KREQIN(:)

INTEGER(KIND=JPIM) :: IPTR,JLEV
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:SETUP_SFLP3',0,zhook_handle)
iptr = ydsc%IPTR
IF(iptr > ydsc%NUMFLDS) THEN
  WRITE(nulerr,*) 'SURFACE FIELDS UNDER-DIMENSINED - GROUP ',&
   & ydsc%CGRPNAME,ydsc%NUMFLDS,kgrib(1),cdname(1)
  CALL abor1('IPTR > YDSC%NUMFLDS')
ENDIF
IF(PRESENT(kgrib)) THEN
  yd%IGRBCODE(:) = kgrib(:)
ENDIF
IF(PRESENT(kreqin)) THEN
  yd%NREQIN(:) = kreqin(:)
ENDIF
IF(PRESENT(cdname)) THEN
  yd%CNAME(:)    = cdname(:)
ENDIF
IF(PRESENT(pdefault)) THEN
  yd%REFVALI(:) = pdefault
ENDIF
IF(PRESENT(ktraj)) THEN
  IF(ktraj == 1) THEN
    DO jlev=1,ydsc%NLEVS
      nstrajgrib(nofftraj+jlev) = yd%IGRBCODE(jlev)
    ENDDO
    nofftraj = nofftraj+ydsc%NLEVS
  ELSEIF(ktraj == 2) THEN
    nofftraj_cst = nofftraj_cst+ydsc%NLEVS
  ELSEIF(ktraj /= 0) THEN
    CALL abor1('SURFACE_FIELDS:SETUP_SFLP3 - UNKNOWN KTRAJ')
  ENDIF
  yd%ITRAJ = ktraj
  ydsc%NDIM5 = ydsc%NDIM5+1
  yd%MP5 = ydsc%NDIM5
ENDIF
DO jlev=1,ydsc%NLEVS
  IF(ydsc%LMTL) THEN
    WRITE(nulout,'(1X,A,2I4,1X,A,6I4)') &
     & ydsc%CGRPNAME(1:6),ydsc%IPTR,jlev,yd%CNAME(jlev),yd%IGRBCODE(jlev),&
     & yd%MP0,yd%MP9,yd%MP1,yd%ITRAJ,yd%NREQIN(jlev)
  ELSE
    WRITE(nulout,'(1X,A,2I4,1X,A,4I4)') &
     & ydsc%CGRPNAME(1:6),ydsc%IPTR,jlev,yd%CNAME(jlev),yd%IGRBCODE(jlev),&
     & yd%MP,yd%ITRAJ,yd%NREQIN(jlev)
  ENDIF
ENDDO
ydsc%IPTR = ydsc%IPTR+1
  
IF (lhook) CALL dr_hook('SURFACE_FIELDS:SETUP_SFLP3',1,zhook_handle)
END SUBROUTINE setup_sflp3

!=========================================================================

SUBROUTINE ini_sflp2(YDSC,YD,KFLDS,LDMTL,CDGRPNAME)
! Initialize 2-D surface field group
TYPE(type_surf_gen),INTENT(INOUT)    :: YDSC
TYPE(type_surf_mtl_2d),INTENT(INOUT) :: YD(:)
INTEGER(KIND=JPIM),INTENT(IN)        :: KFLDS
LOGICAL,INTENT(IN)                   :: LDMTL
CHARACTER(LEN=*),INTENT(IN)          :: CDGRPNAME

INTEGER(KIND=JPIM) :: JFLD, IMAXF
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:INI_SFLP2',0,zhook_handle)

imaxf = SIZE(yd)
ydsc%NUMFLDS = kflds
ydsc%NLEVS = -1
ydsc%IPTR  = 1
ydsc%LMTL  = ldmtl
ydsc%CGRPNAME = cdgrpname
ydsc%NDIM5 = 0
ydsc%NOFFTRAJ = nofftraj 
ydsc%NOFFTRAJ_CST = nofftraj_cst 

nsurf = nsurf+ydsc%NUMFLDS
nsurfl = nsurfl+ydsc%NUMFLDS
IF(ldmtl) THEN
  nprogsurf = nprogsurf+ydsc%NUMFLDS
  nprogsurfl = nprogsurfl+ydsc%NUMFLDS
ENDIF

IF(ldmtl) THEN
  IF (ltwotl) THEN
    ydsc%NDIM = 2*ydsc%NUMFLDS
  ELSE
    ydsc%NDIM = 3*ydsc%NUMFLDS
  ENDIF
ELSE
  ydsc%NDIM = ydsc%NUMFLDS
ENDIF  
ndimsurf = ndimsurf + ydsc%NDIM    
ndimsurfl = ndimsurfl + ydsc%NDIM    
DO jfld=1,kflds
  yd(jfld)%IGRBCODE = -999
  yd(jfld)%CNAME = ''
  yd(jfld)%REFVALI = 0.0_jprb
  yd(jfld)%NREQIN = -1
  yd(jfld)%MP = jfld
  IF (ydsc%LMTL) THEN
    yd(jfld)%MP0 = yd(jfld)%MP
    IF(ltwotl) THEN
      yd(jfld)%MP9 = yd(jfld)%MP0
      yd(jfld)%MP1 = yd(jfld)%MP0+ydsc%NUMFLDS
    ELSE
      yd(jfld)%MP9 = yd(jfld)%MP0+ydsc%NUMFLDS
      yd(jfld)%MP1 = yd(jfld)%MP0+2*ydsc%NUMFLDS
    ENDIF
  ELSE
    yd(jfld)%MP0 = nundefld
    yd(jfld)%MP9 = nundefld
    yd(jfld)%MP1 = nundefld
  ENDIF
  yd(jfld)%MP5 = nundefld
  yd(jfld)%ITRAJ = 0
ENDDO
  
DO jfld=kflds+1,imaxf
  yd(jfld)%MP  = nundefld
  yd(jfld)%MP0 = nundefld
  yd(jfld)%MP9 = nundefld
  yd(jfld)%MP1 = nundefld
  yd(jfld)%MP5 = nundefld
  yd(jfld)%ITRAJ = 0
ENDDO

WRITE(nulout,*) 'INITIALIZING 2-D SURFACE FIELD GROUP ', ydsc%CGRPNAME
WRITE(nulout,*) 'NUMFLDS=',ydsc%NUMFLDS,' LMTL=',ydsc%LMTL 

IF (lhook) CALL dr_hook('SURFACE_FIELDS:INI_SFLP2',1,zhook_handle)
END SUBROUTINE ini_sflp2

!=========================================================================

SUBROUTINE setup_sflp2(YDSC,YD,KGRIB,CDNAME,PDEFAULT,KTRAJ,KREQIN)
! Setup 2-D surface field
TYPE(type_surf_gen),INTENT(INOUT)      :: YDSC
TYPE(type_surf_mtl_2d),INTENT(INOUT)   :: YD
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KGRIB
CHARACTER(LEN=16) ,OPTIONAL,INTENT(IN) :: CDNAME
REAL(KIND=JPRB)   ,OPTIONAL,INTENT(IN) :: PDEFAULT
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KTRAJ
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KREQIN

INTEGER(KIND=JPIM) :: IPTR
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:SETUP_SFLP2',0,zhook_handle)
iptr = ydsc%IPTR
IF(iptr > ydsc%NUMFLDS) THEN
  WRITE(nulerr,*) 'SURFACE FIELDS UNDER-DIMENSINED - GROUP ',ydsc%CGRPNAME,ydsc%NUMFLDS,kgrib,cdname
  CALL abor1('IPTR > YDSC%NUMFLDS')
ENDIF
IF(PRESENT(kgrib)) THEN
  yd%IGRBCODE = kgrib
ENDIF
IF(PRESENT(kreqin)) THEN
  yd%NREQIN = kreqin
ENDIF
IF(PRESENT(cdname)) THEN
  yd%CNAME    = cdname
ENDIF
IF(PRESENT(pdefault)) THEN
  yd%REFVALI = pdefault
ENDIF
IF(PRESENT(ktraj)) THEN
  IF(ktraj == 1) THEN
    nstrajgrib(nofftraj+1) = yd%IGRBCODE
    nofftraj = nofftraj+1
  ELSEIF(ktraj == 2) THEN
    nofftraj_cst = nofftraj_cst+1
  ELSEIF(ktraj /= 0) THEN
    CALL abor1('SURFACE_FIELDS:SETUP_SFLP2 - UNKNOWN KTRAJ')
  ENDIF
  yd%ITRAJ = ktraj
  ydsc%NDIM5 = ydsc%NDIM5+1
  yd%MP5 = ydsc%NDIM5
ENDIF
IF(ydsc%LMTL) THEN
  WRITE(nulout,'(1X,A,I4,1X,A,6I4)') &
   & ydsc%CGRPNAME(1:6),ydsc%IPTR,yd%CNAME,yd%IGRBCODE,&
   & yd%MP0,yd%MP9,yd%MP1,yd%ITRAJ,yd%NREQIN
ELSE
  WRITE(nulout,'(1X,A,I4,1X,A,4I4)') &
   & ydsc%CGRPNAME(1:6),ydsc%IPTR,yd%CNAME,yd%IGRBCODE,yd%MP,yd%ITRAJ,yd%NREQIN
ENDIF
  
ydsc%IPTR = ydsc%IPTR+1
IF (lhook) CALL dr_hook('SURFACE_FIELDS:SETUP_SFLP2',1,zhook_handle)
END SUBROUTINE setup_sflp2

!=========================================================================

SUBROUTINE gppoper(CDACT,KBL,PSP_SB,PSP_SG,PSP_RR,PSP_EP,PSP_X2,YDCOM)
! Operations on prognostic surface fields
CHARACTER(LEN=*),INTENT(IN)            :: CDACT
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KBL
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_SB(:,:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_SG(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_RR(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_EP(:,:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_X2(:,:)
TYPE(type_sfl_comm),OPTIONAL,INTENT(INOUT)   :: YDCOM


REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPPOPER',0,zhook_handle)
IF(PRESENT(kbl)) THEN
  CALL gpoper_3(cdact,sp_sb(:,:,:,kbl),ysp_sbd,ysp_sb%YSB,ydcom)
  CALL gpoper_2(cdact,sp_sg(:,:,kbl)  ,ysp_sgd,ysp_sg%YSG,ydcom)
  CALL gpoper_2(cdact,sp_rr(:,:,kbl)  ,ysp_rrd,ysp_rr%YRR,ydcom)
  CALL gpoper_3(cdact,sp_ep(:,:,:,kbl),ysp_epd,ysp_ep%YEP,ydcom)
  CALL gpoper_2(cdact,sp_x2(:,:,kbl)  ,ysp_x2d,ysp_x2%YX2,ydcom)
ELSE
  IF(PRESENT(psp_sb)) CALL gpoper_3(cdact,psp_sb(:,:,:),ysp_sbd,ysp_sb%YSB,ydcom)
  IF(PRESENT(psp_sg)) CALL gpoper_2(cdact,psp_sg(:,:)  ,ysp_sgd,ysp_sg%YSG,ydcom)
  IF(PRESENT(psp_rr)) CALL gpoper_2(cdact,psp_rr(:,:)  ,ysp_rrd,ysp_rr%YRR,ydcom)
  IF(PRESENT(psp_ep)) CALL gpoper_3(cdact,psp_ep(:,:,:),ysp_epd,ysp_ep%YEP,ydcom)
  IF(PRESENT(psp_x2)) CALL gpoper_2(cdact,psp_x2(:,:)  ,ysp_x2d,ysp_x2%YX2,ydcom)
ENDIF
IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPPOPER',1,zhook_handle)
END SUBROUTINE gppoper

!=========================================================================

SUBROUTINE gpoper(CDACT,KBL,PSP_SB,PSP_SG,PSP_RR,PSD_VF,PSD_VV,YDCOM,PFIELD,PFIELD2)
!Operations on ALL surface groups
CHARACTER(LEN=*),INTENT(IN)            :: CDACT
INTEGER(KIND=JPIM),OPTIONAL,INTENT(IN) :: KBL
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_SB(:,:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_SG(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSP_RR(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSD_VF(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PSD_VV(:,:)
TYPE(type_sfl_comm),OPTIONAL,INTENT(INOUT)   :: YDCOM
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PFIELD(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT) :: PFIELD2(:,:)
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER',0,zhook_handle)
IF(cdact == 'PUTALLFLDS' .OR. cdact == 'GETALLFLDS'   .OR.&
 & cdact == 'TRAJSTORE'  .OR. cdact == 'TRAJSTORECST' .OR. &
 & cdact == 'SET0TOTRAJ' .OR. cdact == 'GETTRAJ'          ) THEN
  IF(.NOT.PRESENT(pfield)) CALL abor1('SURFACE_FIELDS:GPOPER - PFIELD MISSING')
  IF(SIZE(pfield,1) < nproma)  CALL abor1('SURFACE_FIELDS:GPOPER - SIZE(PFIELD,1) < NPROMA)')
ENDIF
IF(cdact == 'PUTALLFLDS' .OR. cdact == 'GETALLFLDS') THEN
  IF(SIZE(pfield,2) < nprogsurfl) CALL abor1('SURFACE_FIELDS:GPOPER - SIZE(PFIELD,2) < NPROGSURFL)')
ENDIF
IF(cdact == 'GETTRAJ') THEN
  IF(.NOT.PRESENT(pfield2)) CALL abor1('SURFACE_FIELDS:GPOPER - PFIELD2 MISSING')
  IF(SIZE(pfield2,1) < nproma)  CALL abor1('SURFACE_FIELDS:GPOPER - SIZE(PFIELD2,1) < NPROMA)')
ENDIF
IF(PRESENT(ydcom)) THEN
  ydcom%L_OK = .false.
  ydcom%IPTRSURF = 0
  ydcom%ICOUNT = 0
ENDIF

nptrsurf = 0
IF(PRESENT(kbl)) THEN
  IF(ysp_sbd%NDIM > 0) THEN
    CALL gpoper_3(cdact,sp_sb(:,:,:,kbl),ysp_sbd,ysp_sb%YSB,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_sgd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sp_sg(:,:,kbl)  ,ysp_sgd,ysp_sg%YSG,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_rrd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sp_rr(:,:,kbl)  ,ysp_rrd,ysp_rr%YRR,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_epd%NDIM > 0) THEN
    CALL gpoper_3(cdact,sp_ep(:,:,:,kbl),ysp_epd,ysp_ep%YEP,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_x2d%NDIM > 0) THEN
    CALL gpoper_2(cdact,sp_x2(:,:,kbl)  ,ysp_x2d,ysp_x2%YX2,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vfd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vf(:,:,kbl)  ,ysd_vfd,ysd_vf%YVF,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vpd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vp(:,:,kbl)  ,ysd_vpd,ysd_vp%YVP,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vvd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vv(:,:,kbl)  ,ysd_vvd,ysd_vv%YVV,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vnd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vn(:,:,kbl)  ,ysd_vnd,ysd_vn%YVN,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vhd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vh(:,:,kbl)  ,ysd_vhd,ysd_vh%YVH,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vad%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_va(:,:,kbl)  ,ysd_vad,ysd_va%YVA,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vcd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vc(:,:,kbl)  ,ysd_vcd,ysd_vc%YVC,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vdd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vd(:,:,kbl)  ,ysd_vdd,ysd_vd%YVD,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_wsd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_ws(:,:,kbl)  ,ysd_wsd,ysd_ws%YWS,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_xad%NDIM > 0) THEN
    CALL gpoper_3(cdact,sd_xa(:,:,:,kbl),ysd_xad,ysd_xa%YXA,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_x2d%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_x2(:,:,kbl)  ,ysd_x2d,ysd_x2%YX2,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vxd%NDIM > 0) THEN
    CALL gpoper_2(cdact,sd_vx(:,:,kbl)  ,ysd_vxd,ysd_vx%YVX,ydcom,pfield,pfield2)
  ENDIF
ELSE
  IF(ysp_sbd%NDIM > 0) THEN
    IF(PRESENT(psp_sb)) &
     & CALL gpoper_3(cdact,psp_sb,ysp_sbd,ysp_sb%YSB,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_sgd%NDIM > 0) THEN
    IF(PRESENT(psp_sg)) &
     & CALL gpoper_2(cdact,psp_sg,ysp_sgd,ysp_sg%YSG,ydcom,pfield,pfield2)
  ENDIF
  IF(ysp_rrd%NDIM > 0) THEN
    IF(PRESENT(psp_rr)) &
     & CALL gpoper_2(cdact,psp_rr,ysp_rrd,ysp_rr%YRR,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vfd%NDIM > 0) THEN
    IF(PRESENT(psd_vf)) &
     & CALL gpoper_2(cdact,psd_vf,ysd_vfd,ysd_vf%YVF,ydcom,pfield,pfield2)
  ENDIF
  IF(ysd_vvd%NDIM > 0) THEN
    IF(PRESENT(psd_vv)) &
     & CALL gpoper_2(cdact,psd_vv,ysd_vvd,ysd_vv%YVV,ydcom,pfield,pfield2)
  ENDIF
ENDIF
IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER',1,zhook_handle)
END SUBROUTINE gpoper

!=========================================================================

SUBROUTINE gpoper_2(CDACT,PFLD,YDSC,YD,YDCOM,PFIELD,PFIELD2)
! Operations on 2-D surface groups
CHARACTER(LEN=*),INTENT(IN)                :: CDACT
REAL(KIND=JPRB),INTENT(INOUT)              :: PFLD(:,:)
TYPE(type_surf_gen),INTENT(IN)             :: YDSC
TYPE(type_surf_mtl_2d),INTENT(IN)          :: YD(:)
TYPE(type_sfl_comm),OPTIONAL,INTENT(INOUT) :: YDCOM
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT)     :: PFIELD(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT)     :: PFIELD2(:,:)

INTEGER(KIND=JPIM) :: J,IPTR,IPTR2
REAL(KIND=JPRB) :: ZZPHY
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER_2',0,zhook_handle)
IF(cdact == 'SET9TO0') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP9) = pfld(:,yd(j)%MP0)
  ENDDO
ELSEIF(cdact == 'SET1TO0') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP1) = pfld(:,yd(j)%MP0)
  ENDDO
ELSEIF(cdact == 'SET1TO9') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP1) = pfld(:,yd(j)%MP9)
  ENDDO
ELSEIF(cdact == 'SET1TO9AD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP0) = pfld(:,yd(j)%MP9)+pfld(:,yd(j)%MP1)
    pfld(:,yd(j)%MP1) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET0TO1') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP0) = pfld(:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'SET0TO1AD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP1) = pfld(:,yd(j)%MP1)+pfld(:,yd(j)%MP0)
    pfld(:,yd(j)%MP0) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET9TO1') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP9) = pfld(:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'PHTFILT') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  zzphy=1.0_jprb-repsp1
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP9) = repsp1*pfld(:,yd(j)%MP1)+zzphy*pfld(:,yd(j)%MP0)
    pfld(:,yd(j)%MP0) = pfld(:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'PHTFILTAD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  zzphy=1.0_jprb-repsp1
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP1) = pfld(:,yd(j)%MP1)+pfld(:,yd(j)%MP0)
    pfld(:,yd(j)%MP0) = 0.0_jprb
    pfld(:,yd(j)%MP1) = pfld(:,yd(j)%MP1)+repsp1*pfld(:,yd(j)%MP9)
    pfld(:,yd(j)%MP0) = pfld(:,yd(j)%MP0)+zzphy *pfld(:,yd(j)%MP9)
    pfld(:,yd(j)%MP9) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET0TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP0) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SET9TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP9) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SET1TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_2 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,yd(j)%MP1) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SETALLTOVAL') THEN
  DO j=1,ydsc%NDIM
    pfld(:,j) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SETDEFAULT') THEN
  DO j=1,ydsc%NUMFLDS
    IF(yd(j)%NREQIN == -1) THEN
      pfld(:,yd(j)%MP) = yd(j)%REFVALI
    ENDIF
  ENDDO
ELSEIF(cdact == 'TRAJSTORE') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        iptr = iptr+1
        pfield(:,iptr) = pfld(:,yd(j)%MP)
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'TRAJSTORECST') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr2 = ydsc%NOFFTRAJ_CST
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 2) THEN
        iptr2 = iptr2+1
        pfield(:,iptr2) = pfld(:,yd(j)%MP)
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'SET0TOTRAJ') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        iptr = iptr+1
        pfld(:,yd(j)%MP) = pfield(:,iptr)
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'GETTRAJ') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    iptr2 = ydsc%NOFFTRAJ_CST
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        iptr = iptr+1
        pfld(:,yd(j)%MP5) = pfield(:,iptr)
      ELSEIF(yd(j)%ITRAJ == 2) THEN
        iptr2 = iptr2+1
        pfld(:,yd(j)%MP5) = pfield2(:,iptr2)
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'GETALLFLDS') THEN
  DO j=1,ydsc%NDIM
    nptrsurf = nptrsurf+1
    pfield(:,nptrsurf) = pfld(:,j)
  ENDDO
ELSEIF(cdact == 'PUTALLFLDS') THEN
  DO j=1,ydsc%NDIM
    nptrsurf = nptrsurf+1
    pfld(:,j) = pfield(:,nptrsurf)
  ENDDO
ELSEIF(cdact == 'GETGRIBPOS') THEN
  DO j=1,ydsc%NUMFLDS  
    ydcom%IPTRSURF = ydcom%IPTRSURF+1
    IF(yd(j)%IGRBCODE == ydcom%IGRBCODE) THEN
      ydcom%IFLDNUM  = ydcom%IPTRSURF
      ydcom%L_OK = .true.
    ENDIF
  ENDDO
ELSEIF(cdact == 'GETFIELD') THEN
  DO j=1,ydsc%NUMFLDS  
    ydcom%IPTRSURF = ydcom%IPTRSURF+1
    IF(ydcom%IPTRSURF == ydcom%IFLDNUM) THEN
      pfield(:,1) = pfld(:,j)
      ydcom%L_OK = .true.
    ENDIF
  ENDDO
ELSEIF(cdact == 'GRIBIN') THEN
  DO j=1,ydsc%NUMFLDS  
    ydcom%IPTRSURF = ydcom%IPTRSURF+1
    IF(yd(j)%NREQIN == 1) THEN
      ydcom%ICOUNT = ydcom%ICOUNT+1
      ydcom%ICODES(ydcom%ICOUNT) = yd(j)%IGRBCODE
    ENDIF
  ENDDO
ELSE
  WRITE(nulout,*) 'SURFACE_FIELD:GPPOPER UNKNOWN ACTION - ',cdact
  CALL abor1('SURFACE_FIELD:GPPOPER - UNKNOWN ACTION')
ENDIF
IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER_2',1,zhook_handle)
END SUBROUTINE gpoper_2

!=========================================================================

SUBROUTINE gpoper_3(CDACT,PFLD,YDSC,YD,YDCOM,PFIELD,PFIELD2)
! Operations on 3-D surface groups
CHARACTER(LEN=*),INTENT(IN)                :: CDACT
REAL(KIND=JPRB),INTENT(INOUT)              :: PFLD(:,:,:)
TYPE(type_surf_gen),INTENT(IN)             :: YDSC
TYPE(type_surf_mtl_3d),INTENT(IN)          :: YD(:)
TYPE(type_sfl_comm),OPTIONAL,INTENT(INOUT) :: YDCOM
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT)     :: PFIELD(:,:)
REAL(KIND=JPRB),OPTIONAL,INTENT(INOUT)     :: PFIELD2(:,:)

INTEGER(KIND=JPIM) :: J,JLEV,IPTR,IPTR2
REAL(KIND=JPRB) :: ZZPHY
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!-------------------------------------------------------------------------

IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER_3',0,zhook_handle)
IF(cdact == 'SET9TO0') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP9) = pfld(:,:,yd(j)%MP0)
  ENDDO
ELSEIF(cdact == 'SET1TO0') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP1) = pfld(:,:,yd(j)%MP0)
  ENDDO
ELSEIF(cdact == 'SET1TO9') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP1) = pfld(:,:,yd(j)%MP9)
  ENDDO
ELSEIF(cdact == 'SET1TO9AD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP9) = pfld(:,:,yd(j)%MP9)+pfld(:,:,yd(j)%MP1)
    pfld(:,:,yd(j)%MP1) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET0TO1') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP0) = pfld(:,:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'SET0TO1AD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP1) = pfld(:,:,yd(j)%MP1)+pfld(:,:,yd(j)%MP0)
    pfld(:,:,yd(j)%MP0) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET9TO1') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP9) = pfld(:,:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'PHTFILT') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  zzphy=1.0_jprb-repsp1
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP9) = repsp1*pfld(:,:,yd(j)%MP1)+zzphy*pfld(:,:,yd(j)%MP0)
    pfld(:,:,yd(j)%MP0) = pfld(:,:,yd(j)%MP1)
  ENDDO
ELSEIF(cdact == 'PHTFILTAD') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  zzphy=1.0_jprb-repsp1
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP1) = pfld(:,:,yd(j)%MP1)+pfld(:,:,yd(j)%MP0)
    pfld(:,:,yd(j)%MP0) = 0.0_jprb
    pfld(:,:,yd(j)%MP1) = pfld(:,:,yd(j)%MP1)+repsp1*pfld(:,:,yd(j)%MP9)
    pfld(:,:,yd(j)%MP0) = pfld(:,:,yd(j)%MP0)+zzphy *pfld(:,:,yd(j)%MP9)
    pfld(:,:,yd(j)%MP9) = 0.0_jprb
  ENDDO
ELSEIF(cdact == 'SET0TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP0) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SET9TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP9) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SET1TOVAL') THEN
  IF( .NOT. ydsc%LMTL) CALL abor1('SURFACE_FIELDS:GPOPER_3 : FIELD NOT MULTI-TIME LEVEL')
  DO j=1,ydsc%NUMFLDS
    pfld(:,:,yd(j)%MP1) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SETALLTOVAL') THEN
  DO j=1,ydsc%NDIM
    pfld(:,:,j) = ydcom%VALUE
  ENDDO
ELSEIF(cdact == 'SETDEFAULT') THEN
  DO j=1,ydsc%NUMFLDS
    DO jlev=1,ydsc%NLEVS
      IF(yd(j)%NREQIN(jlev) == -1) THEN
        pfld(:,jlev,yd(j)%MP) = yd(j)%REFVALI(jlev)
      ENDIF
    ENDDO
  ENDDO
ELSEIF(cdact == 'TRAJSTORE') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        DO jlev=1,ydsc%NLEVS
          iptr = iptr+1
          pfield(:,iptr) = pfld(:,jlev,yd(j)%MP)
        ENDDO
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'TRAJSTORECST') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr2 = ydsc%NOFFTRAJ_CST
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 2) THEN
        DO jlev=1,ydsc%NLEVS
          iptr2 = iptr2+1
          pfield(:,iptr2) = pfld(:,jlev,yd(j)%MP)
        ENDDO
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'SET0TOTRAJ') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        DO jlev=1,ydsc%NLEVS
          iptr = iptr+1
          pfld(:,jlev,yd(j)%MP) = pfield(:,iptr)
        ENDDO
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'GETTRAJ') THEN
  IF(ydsc%NDIM5 > 0 ) THEN
    iptr = ydsc%NOFFTRAJ
    iptr2 = ydsc%NOFFTRAJ_CST
    DO j=1,ydsc%NUMFLDS
      IF(yd(j)%ITRAJ == 1) THEN
        DO jlev=1,ydsc%NLEVS
          iptr = iptr+1
          pfld(:,jlev,yd(j)%MP5) = pfield(:,iptr)
        ENDDO
      ELSEIF(yd(j)%ITRAJ == 2) THEN
        DO jlev=1,ydsc%NLEVS
          iptr2 = iptr2+1
          pfld(:,jlev,yd(j)%MP5) = pfield2(:,iptr2)
        ENDDO
      ENDIF
    ENDDO
  ENDIF
ELSEIF(cdact == 'GETALLFLDS') THEN
  DO j=1,ydsc%NDIM
    DO jlev=1,ydsc%NLEVS
      nptrsurf = nptrsurf+1
      pfield(:,nptrsurf) = pfld(:,jlev,j)
    ENDDO
  ENDDO
ELSEIF(cdact == 'PUTALLFLDS') THEN
  DO j=1,ydsc%NDIM
    DO jlev=1,ydsc%NLEVS
      nptrsurf = nptrsurf+1
      pfld(:,jlev,j) = pfield(:,nptrsurf)
    ENDDO
  ENDDO
ELSEIF(cdact == 'GETGRIBPOS') THEN
  DO j=1,ydsc%NUMFLDS  
    DO jlev=1,ydsc%NLEVS
      ydcom%IPTRSURF = ydcom%IPTRSURF+1
      IF(yd(j)%IGRBCODE(jlev) == ydcom%IGRBCODE) THEN
        ydcom%IFLDNUM  = ydcom%IPTRSURF
        ydcom%L_OK = .true.
      ENDIF
    ENDDO
  ENDDO
ELSEIF(cdact == 'GETFIELD') THEN
  DO j=1,ydsc%NUMFLDS  
    DO jlev=1,ydsc%NLEVS
      ydcom%IPTRSURF = ydcom%IPTRSURF+1
      IF(ydcom%IPTRSURF == ydcom%IFLDNUM) THEN
        pfield(:,1) = pfld(:,jlev,j)
        ydcom%L_OK = .true.
      ENDIF
    ENDDO
  ENDDO
ELSEIF(cdact == 'GRIBIN') THEN
  DO j=1,ydsc%NUMFLDS  
    DO jlev=1,ydsc%NLEVS
      ydcom%IPTRSURF = ydcom%IPTRSURF+1
      IF(yd(j)%NREQIN(jlev) == 1) THEN
        ydcom%ICOUNT = ydcom%ICOUNT+1
        ydcom%ICODES(ydcom%ICOUNT) = yd(j)%IGRBCODE(jlev)
      ENDIF
    ENDDO
  ENDDO
ELSE
  WRITE(nulout,*) 'SURFACE_FIELD:GPPOPER UNKNOWN ACTION - ',cdact
  CALL abor1('SURFACE_FIELD:GPPOPER - UNKNOWN ACTION')
ENDIF
IF (lhook) CALL dr_hook('SURFACE_FIELDS:GPOPER_3',1,zhook_handle)
END SUBROUTINE gpoper_3

!=========================================================================

SUBROUTINE surf_store
! Store all surface fields
INTEGER(KIND=JPIM) :: JBL
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('SURFACE_FIELDS:SURF_STORE',0,zhook_handle)
ALLOCATE(surf_store_array(nproma,ndimsurfl,ngpblks))
DO jbl=1,ngpblks
  CALL gpoper('GETALLFLDS',kbl=jbl,pfield=surf_store_array(:,:,jbl))
ENDDO
IF (lhook) CALL dr_hook('SURFACE_FIELDS:SURF_STORE',1,zhook_handle)
END SUBROUTINE surf_store

!=========================================================================

SUBROUTINE surf_restore
! Restore all surface fields
INTEGER(KIND=JPIM) :: JBL
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('SURFACE_FIELDS:SURF_RESTORE',0,zhook_handle)
IF(.NOT. ALLOCATED(surf_store_array)) &
 & CALL abor1('SURFACE_FIELDS:SURF_RESTORE - SURF_STORE NOT ALLOCATED')
DO jbl=1,ngpblks
  CALL gpoper('PUTALLFLDS',kbl=jbl,pfield=surf_store_array(:,:,jbl))
ENDDO
DEALLOCATE(surf_store_array)
IF (lhook) CALL dr_hook('SURFACE_FIELDS:SURF_RESTORE',1,zhook_handle)

END SUBROUTINE surf_restore

!=========================================================================

SUBROUTINE allo_surf
! Allocate surface field arrays
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('SURFACE_FIELDS:ALLO_SURF',0,zhook_handle)
ALLOCATE(sp_sb(nproma,ysp_sbd%NLEVS,ysp_sbd%NDIM,ngpblks))
ALLOCATE(sp_sg(nproma,ysp_sgd%NDIM,ngpblks))
ALLOCATE(sp_rr(nproma,ysp_rrd%NDIM,ngpblks))
ALLOCATE(sp_ep(nproma,ysp_epd%NLEVS,ysp_epd%NDIM,ngpblks))
ALLOCATE(sp_x2(nproma,ysp_x2d%NDIM,ngpblks))
ALLOCATE(sp_ci(nproma,ysp_cid%NDIM,ngpblks))
ALLOCATE(sd_vf(nproma,ysd_vfd%NDIM,ngpblks))
ALLOCATE(sd_vp(nproma,ysd_vpd%NDIM,ngpblks))
ALLOCATE(sd_vv(nproma,ysd_vvd%NDIM,ngpblks))
ALLOCATE(sd_vn(nproma,ysd_vnd%NDIM,ngpblks))
ALLOCATE(sd_vh(nproma,ysd_vhd%NDIM,ngpblks))
ALLOCATE(sd_va(nproma,ysd_vad%NDIM,ngpblks))
ALLOCATE(sd_vc(nproma,ysd_vcd%NDIM,ngpblks))
ALLOCATE(sd_vd(nproma,ysd_vdd%NDIM,ngpblks))
ALLOCATE(sd_ws(nproma,ysd_wsd%NDIM,ngpblks))
ALLOCATE(sd_xa(nproma,ysd_xad%NLEVS,ysd_xad%NDIM,ngpblks))
ALLOCATE(sd_x2(nproma,ysd_x2d%NDIM,ngpblks))
ALLOCATE(sd_vx(nproma,ysd_vxd%NDIM,ngpblks))

IF (lhook) CALL dr_hook('SURFACE_FIELDS:ALLO_SURF',1,zhook_handle)
END SUBROUTINE allo_surf

!=========================================================================

SUBROUTINE deallo_surf
! Deallocate surface field arrays
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('SURFACE_FIELDS:DEALLO_SURF',0,zhook_handle)
IF(ALLOCATED(sp_sb)) DEALLOCATE(sp_sb)
IF(ALLOCATED(sp_sg)) DEALLOCATE(sp_sg)
IF(ALLOCATED(sp_rr)) DEALLOCATE(sp_rr)
IF(ALLOCATED(sp_ep)) DEALLOCATE(sp_ep)
IF(ALLOCATED(sp_x2)) DEALLOCATE(sp_x2)
IF(ALLOCATED(sp_ci)) DEALLOCATE(sp_ci)
IF(ALLOCATED(sd_vf)) DEALLOCATE(sd_vf)
IF(ALLOCATED(sd_vp)) DEALLOCATE(sd_vp)
IF(ALLOCATED(sd_vv)) DEALLOCATE(sd_vv)
IF(ALLOCATED(sd_vn)) DEALLOCATE(sd_vn)
IF(ALLOCATED(sd_vh)) DEALLOCATE(sd_vh)
IF(ALLOCATED(sd_va)) DEALLOCATE(sd_va)
IF(ALLOCATED(sd_vc)) DEALLOCATE(sd_vc)
IF(ALLOCATED(sd_vd)) DEALLOCATE(sd_vd)
IF(ALLOCATED(sd_ws)) DEALLOCATE(sd_ws)
IF(ALLOCATED(sd_xa)) DEALLOCATE(sd_xa)
IF(ALLOCATED(sd_x2)) DEALLOCATE(sd_x2)
IF(ALLOCATED(sd_vx)) DEALLOCATE(sd_vx)
IF (lhook) CALL dr_hook('SURFACE_FIELDS:DEALLO_SURF',1,zhook_handle)
END SUBROUTINE deallo_surf

!=========================================================================

END MODULE surface_fields