ol_time_interp_atm.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.
!     ######spl
SUBROUTINE ol_time_interp_atm (KSURF_STEP,KNB_ATM,             &
                               PTA1,PTA2,PQA1,PQA2,PWIND1,PWIND2,        &
                               PDIR_SW1,PDIR_SW2,PSCA_SW1,PSCA_SW2,      &
                               PLW1,PLW2,PSNOW2,PRAIN2,                  &
                               PPS1,PPS2,PCO21,PCO22,PDIR1,PDIR2,        &
                               PZEN,PSUMZEN )  
!**************************************************************************
!
!!    PURPOSE
!!    -------
!        Time interpolation of the atmospheric forcing
!        So far, it is a simple linear interpolation.
!        More complex interpolation may be added, especially for the atmospheric
!        radiation (Option to use).
!        Output are in the module 
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      F. Habets   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/2003
!
USE modn_io_offline, ONLY : linterp_sw
!
USE modd_csts,       ONLY : xpi, xrd, xrv, xg
USE modd_surf_par,   ONLY : xundef
USE modd_forc_atm,  ONLY: xta         ,&! air temperature forcing               (K)
                            xqa       ,&! air specific humidity forcing         (kg/m3)
                            xrhoa     ,&! air density forcing                   (kg/m3)
                            xzs       ,&! orography                             (m)
                            xu        ,&! zonal wind                            (m/s)
                            xv        ,&! meridian wind                         (m/s)
                            xdir_sw   ,&! direct  solar radiation (on horizontal surf.)
                            xsca_sw   ,&! diffuse solar radiation (on horizontal surf.)
                            xlw       ,&! longwave radiation (on horizontal surf.)
                            xps       ,&! pressure at atmospheric model surface (Pa)
                            xpa       ,&! pressure at forcing level             (Pa)
                            xrhoa     ,&! density at forcing level              (kg/m3)
                            xco2      ,&! CO2 concentration in the air          (kg/kg)
                            xsnow     ,&! snow precipitation                    (kg/m2/s)
                            xrain     ,&! liquid precipitation                  (kg/m2/s)
                            xzref       ! height of T,q forcing                 (m)  
!
USE modi_get_luout
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64
!$ USE OMP_LIB
#endif
!
IMPLICIT NONE
!
!
#ifndef AIX64
!$ INCLUDE 'omp_lib.h'
#endif
!
! global variables
INTEGER,INTENT(IN) :: KSURF_STEP, KNB_ATM
REAL, DIMENSION(:),INTENT(IN) :: PTA1,PTA2,PQA1,PQA2,PWIND1,PWIND2
REAL, DIMENSION(:),INTENT(IN) :: PDIR_SW1,PDIR_SW2,PSCA_SW1,PSCA_SW2,PLW1,PLW2
REAL, DIMENSION(:),INTENT(IN) :: PSNOW2,PRAIN2,PPS1,PPS2,PCO21,PCO22,PDIR1,PDIR2
REAL, DIMENSION(:),INTENT(IN) :: PZEN,PSUMZEN 

! local variables
REAL :: ZDTA, ZDQA, ZDDIR_SW, ZDSCA_SW, ZDLW,  &
        ZDPS, ZDCO2, ZDU, ZDV, ZU1, ZV1, ZU2, ZV2 
REAL :: ZPI, ZNB_ATM, ZSURF_STEP, ZCOEF, ZCOEF2
INTEGER :: J
INTEGER :: ILUOUT
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!========================================================================
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_1',0,zhook_handle)
!
 CALL get_luout('OFFLIN',iluout)
!
zpi = xpi/180.
znb_atm = knb_atm*1.
zsurf_step = ksurf_step*1.-1.
zcoef = zsurf_step / znb_atm
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_1',1,zhook_handle)
!
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP)
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_2',0,zhook_handle_omp)
!$OMP DO PRIVATE(J,ZU1,ZU2,ZV1,ZV2,ZDU,ZDV,ZDTA, &
!$OMP ZDQA,ZDLW,ZDPS,ZDCO2,ZDDIR_SW,ZDSCA_SW,ZCOEF2)
DO j = 1,SIZE(pta1)
  !
  IF (pta1(j)/=xundef) THEN
    ! 
    ! Compute wind components
    !
    ! zonal wind
    zu1 = pwind1(j) * sin(pdir1(j)*zpi)
    zu2 = pwind2(j) * sin(pdir2(j)*zpi)
    zdu = (zu2-zu1)*zcoef
    xu(j) = zu1 + zdu
    !
    zv1 = pwind1(j) * cos(pdir1(j)*zpi) 
    zv2 = pwind2(j) * cos(pdir2(j)*zpi)
    zdv = (zv2-zv1)*zcoef
    xv(j) = zv1 + zdv
    !    
    ! Compute variation from atmospheric time step J and J+1
    !
    zdta     = (pta2(j)-pta1(j))*zcoef
    xta(j) = pta1(j) + zdta
    !
    zdqa     = (pqa2(j)-pqa1(j))*zcoef
    xqa(j) = pqa1(j) + zdqa
    !
    zdlw     = (plw2(j)-plw1(j))*zcoef
    xlw(j) = plw1(j) + zdlw
    !
    zdps     = (pps2(j)-pps1(j))*zcoef
    xps(j) = pps1(j) + zdps
    !
    zdco2    = (pco22(j)-pco21(j))*zcoef
    xco2(j) = pco21(j) + zdco2
    !
    IF (linterp_sw) THEN
      !
      zcoef2=0.
      IF (psumzen(j)>0.) zcoef2=max((cos(pzen(j))/psumzen(j)),0.)
      !
      xdir_sw(j,1) = min(pdir_sw2(j)*zcoef2,1300.0*max(cos(pzen(j)),0.))
      !
      xsca_sw(j,1) = min(psca_sw2(j)*zcoef2,1300.0*max(cos(pzen(j)),0.))
      !
    ELSE
      !
      zddir_sw = (pdir_sw2(j)-pdir_sw1(j))*zcoef
      xdir_sw(j,1) = pdir_sw1(j)+zddir_sw
      !
      zdsca_sw = (psca_sw2(j)-psca_sw1(j))*zcoef
      xsca_sw(j,1) = psca_sw1(j)+zdsca_sw
      !
    ENDIF
    !
    !
    xrain(j)= prain2(j)
    xsnow(j)= psnow2(j)
    !
    !
    xrhoa(j) = xps(j) / ( xta(j)*xrd * ( 1.+((xrv/xrd)-1.)*xqa(j) ) + xzref(j)*xg )
    !
    ! humidity in kg/m3
    xqa(j) = xqa(j) * xrhoa(j)
    !
  ENDIF
  !
ENDDO
!$OMP END DO
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_2',1,zhook_handle_omp)
!$OMP END PARALLEL
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_3',0,zhook_handle)
! air density
!
! Check No value data
!---------------------
! Error cases
!
IF ((minval(xta)  .EQ.xundef).OR.(minval(xqa).EQ.xundef).OR.&
      (minval(xu).EQ.xundef).OR.(minval(xrain).EQ.xundef).OR.&
      (minval(xsnow).EQ.xundef)) THEN  
    WRITE(iluout,*)'MINVAL(XTA),MINVAL(XQA),MINVAL(XU),MINVAL(XRAIN),MINVAL(XSNOW)'
    WRITE(iluout,*)minval(xta),minval(xqa),minval(xu),minval(xrain),minval(xsnow)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF ((minval(xdir_sw).EQ.xundef).AND.(minval(xsca_sw).EQ.xundef)) THEN
    WRITE(iluout,*)'MINVAL(XSCA_SW),MINVAL(XDIR_SW)'
    WRITE(iluout,*)minval(xsca_sw),minval(xdir_sw)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF ((minval(xps).EQ.xundef).AND.(minval(xzs).EQ.xundef)) THEN
    WRITE(iluout,*)'MINVAL(XPS),MINVAL(XZS)'
    WRITE(iluout,*)minval(xps),minval(xzs)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF (minval(xdir_sw).EQ.xundef) xdir_sw(:,:)=0. ! No direct solar radiation
IF (minval(xsca_sw).EQ.xundef) xsca_sw(:,:)=0. ! No diffuse solar radiation
IF (minval(xps)    .EQ.xundef) THEN            ! No surface Pressure 
   WRITE(iluout,*)' OL_TIME_INTERP_ATM: SURFACE PRESSURE COMPUTED FROM ZS'
   xps(:)  = 101325*(1-0.0065 * xzs(:)/288.15)**5.31
ENDIF
!
!* forcing level pressure from hydrostatism
WHERE(xps(:)/=xundef)
  xpa(:) = xps(:) - xrhoa(:) * xzref(:) * xg
ENDWHERE
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM_3',1,zhook_handle)
!
END SUBROUTINE ol_time_interp_atm