sso.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 sso (UG, USS, &
                      osso,osso_anis,psea)
!     #########################
!
!!*SSO  computes the SSO anisotropy, direction and slope
!!
!!
!!    METHOD
!!    ------
!!    See Lott and Miller, 1997, QJRMS 101-127
!!   
!!    AUTHOR
!!    ------
!!
!!    V.Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    23/05/97
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modi_get_mesh_dim
USE modi_get_adjacent_meshes
!
USE modd_csts,           ONLY : xpi
USE modd_pgdwork,        ONLY : nsso, xssqo, lssqo
USE modd_pgd_grid,       ONLY : nl, cgrid, xgrid_par, ngrid_par
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of dummy arguments
!            ------------------------------
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_sso_t), INTENT(INOUT) :: uss
!
LOGICAL, DIMENSION(:), INTENT(OUT) :: osso   ! .T. : the SSO coefficients
!                                            ! are computed at grid point
!                                            ! .F. : not enough sub-grid
!                                            ! information avalaible to
!                                            ! compute the coefficients
LOGICAL, DIMENSION(:), INTENT(OUT) :: osso_anis ! .T. : the SSO anisotropy
!                                            ! are computed at grid point
!                                            ! .F. : not enough sub-grid
!                                            ! information avalaible to
!                                            ! compute the coefficients
REAL,    DIMENSION(:), INTENT(IN)  :: psea   ! sea fraction
!
!*    0.2    Declaration of indexes
!            ----------------------
!
INTEGER :: jl          ! loop index on grid meshs
INTEGER :: il          ! grid mesh index of second subgrid point used
INTEGER :: jiss, jjss  ! loop indexes for subsquares arrays
INTEGER :: jnext       ! loop index on subgrid meshes
INTEGER :: jprev       ! loop index on subgrid meshes
INTEGER :: inext       ! index to add to JISS or JJSS to obtain the following
!                      ! point containing a data in a segment
INTEGER :: iprev       ! index to remove from JISS or JJSS to obtain the
!                      ! previous point containing a data in a segment
!
INTEGER :: imaxi       ! index of the next subsquare with data along I axis,
                       ! or last subsquare inside grid mesh along I axis
INTEGER :: imaxj       ! index of the next subsquare with data along J axis,
                       ! or last subsquare inside grid mesh along J axis
INTEGER, DIMENSION(NL) :: ileft   ! index of left   grid mesh 
INTEGER, DIMENSION(NL) :: iright  ! index of right  grid mesh 
INTEGER, DIMENSION(NL) :: itop    ! index of top    grid mesh 
INTEGER, DIMENSION(NL) :: ibottom ! index of bottom grid mesh
!
!*    0.3    Declaration of working arrays inside a MESONH grid (JI,JJ)
!            -----------------------------
!
REAL,    DIMENSION(NSSO,NSSO) :: zdzsdx  ! topographic gradient along x
REAL,    DIMENSION(NSSO,NSSO) :: zdzsdy  ! topographic gradient along y
LOGICAL, DIMENSION(NSSO,NSSO) :: gdzsdx  ! occurence of gradient along x
LOGICAL, DIMENSION(NSSO,NSSO) :: gdzsdy  ! occurence of gradient along y
REAL :: zdxeff    ! width of a subsquare along I axis
REAL :: zdyeff    ! width of a subsquare along J axis
LOGICAL :: gbound ! .T.: first left (for dzs/dx) or first bottom (for dzs/dy)
                  ! sub-square gradients is being computed
!
!*    0.4    Declaration of other local variables
!            ------------------------------------
!
REAL, DIMENSION(NL) :: zdx        ! grid mesh size in x direction
REAL, DIMENSION(NL) :: zdy        ! grid mesh size in y direction
REAL, DIMENSION(NL) :: zhxx       ! topographic gradient correlation tensor: x,x
REAL, DIMENSION(NL) :: zhxy       ! topographic gradient correlation tensor: x,y
REAL, DIMENSION(NL) :: zhyy       ! topographic gradient correlation tensor: y,y
REAL, DIMENSION(NL) :: zk, zl, zm ! diagonalised terms of the tensor
REAL(KIND=JPRB) :: zhook_handle
!
!----------------------------------------------------------------------------
!
!*    1.     Initializations
!            ---------------
!
IF (lhook) CALL dr_hook('SSO',0,zhook_handle)
zk=0.
zl=0.
zm=0.
!
!*    1.1    Occurence of computation of the coefficients
!            --------------------------------------------
!
osso(:)=.false.
osso_anis(:)=.false.
!
!
!*    1.2    Grid dimension (meters)
!            -----------------------
!
 CALL get_mesh_dim(cgrid,ngrid_par,nl,xgrid_par,zdx,zdy,ug%XMESH_SIZE)
!
!
!*    1.3    Left, top, right and bottom adjacent gris meshes
!            ------------------------------------------------
!
 CALL get_adjacent_meshes(cgrid,ngrid_par,nl,xgrid_par,ileft,iright,itop,ibottom)
!
!----------------------------------------------------------------------------
!
!*    2.     Loop on MESONH grid points
!            --------------------------
!
zhyy(:) = 0.
zhxx(:) = 0.
zhxy(:) = 0.
!
DO jl=1,nl
!
!
!*    2.1    No land in grid mesh
!            --------------------
!
    IF (psea(jl)==1.) cycle
!
    zdxeff=zdx(jl)/float(nsso)
    zdyeff=zdy(jl)/float(nsso)
!
!*    2.3    Not enough data for computation
!            -------------------------------
!
! 1st step: removes points where orography data is not present at all.
!
    IF ( count( lssqo(:,:,jl)) == 0  ) cycle
!
!----------------------------------------------------------------------------
!
!*    3.     Computations of the gradients along x
!            -------------------------------------
!
    gdzsdx(:,:)=.false.
    zdzsdx(:,:)= 0.
!
!*    3.1    loop on jss index (there is no specific computation along j)
!            -----------------
!
    DO jjss=1,nsso
!
!* left point mark initialization
!
        gbound=.true.
!
!*    3.2    loop on iss index
!            -----------------
!
      DO jiss=1,nsso
!
!
!*    3.3    search for two consecutive grid points
!            --------------------------------------
!
!*    3.3.1 first one
!
        IF (.NOT. lssqo(jiss,jjss,jl) ) cycle
!
!*    3.3.2  second one (up to one grid mesh further)
!
        DO jnext=1,nsso
          IF (jiss+jnext>nsso) THEN
            il    = iright(jl)
            inext = jiss+jnext-nsso
          ELSE
            il    = jl
            inext = jiss+jnext
          END IF
          ! no right point
          IF (il==0) EXIT
          ! subgrid data found
          IF (lssqo(inext,jjss,il)) EXIT
        END DO
!
!*    3.3.3  none found: end of loop along jss
!                        ---------------------
!
        IF (jnext>=nsso+1) EXIT
!
!*    3.3.4  second point outside of the domain: end of loop along iss
!                                                ---------------------
!
        IF (il==0) EXIT
!
!*    3.4    dzs/dx term
!            -----------
!
        imaxi = min(jiss+jnext-1,nsso)
!
        zdzsdx(jiss:imaxi,jjss) = ( xssqo(inext,jjss,il) - xssqo(jiss,jjss,jl)) &
                                    / float(jnext) / zdxeff  
!
        gdzsdx(jiss:imaxi,jjss) = .true.
!
!
!*    3.5    left data point not on the left of the grid mesh (one more computation)
!            ------------------------------------------------
!
        IF (gbound .AND. jiss/=1) THEN
          DO jprev=1,nsso
            IF (jiss-jprev<1) THEN
              il    = ileft(jl)
              iprev = jiss-jprev+nsso
            ELSE
              il    = jl
              iprev = jiss-jprev
            END IF
            ! no left point
             IF (il==0) EXIT
            ! subgrid data found
            IF (lssqo(iprev,jjss,il)) EXIT
          END DO

          IF (.NOT. (jprev>=nsso+1 .OR. il==0)) THEN
            zdzsdx(1:jiss,jjss) = ( xssqo(jiss,jjss,jl) - xssqo(iprev,jjss,il)) &
                                    / float(jprev) / zdxeff  
!
            gdzsdx(1:jiss,jjss) = .true.
          END IF
        END IF
!
!
        gbound=.false.
!
!
!*    3.6    end of loop on iss index
!            ------------------------
!
      END DO
!
!*    3.7    end of loop on jss index
!            ------------------------
!
    END DO
!
!----------------------------------------------------------------------------
!
!*    4.     Computations of the gradients along y
!            -------------------------------------
!
    gdzsdy(:,:)=.false.
    zdzsdy(:,:)= 0.
!
!*    4.1    loop on iss index (there is no specific computation along i)
!            -----------------
!
    DO jiss=1,nsso
!
!* bottom point mark initialization
!
      gbound=.true.
!
!*    4.2    loop on jss index
!            -----------------
!
      DO jjss=1,nsso
!
!
!*    4.3    search for two consecutive grid points
!            --------------------------------------
!
!*    4.3.1 first one
!
        IF (.NOT. lssqo(jiss,jjss,jl) ) cycle
!
!*    4.3.2  second one (up to one grid mesh further)
!
        DO jnext=1,nsso
          IF (jjss+jnext>nsso) THEN
            il    = itop(jl)
            inext = jjss+jnext-nsso
          ELSE
            il    = jl
            inext = jjss+jnext
          END IF
          ! no top point
          IF (il==0) EXIT
          ! subgrid data found
          IF (lssqo(jiss,inext,il)) EXIT
        END DO
!
!*    4.3.3  none found: end of loop along iss
!                        ---------------------
!
        IF (jnext>=nsso+1) EXIT
!
!*    3.3.4  second point outside of the domain: end of loop along iss
!                                                ---------------------
!
        IF (il==0) EXIT
!
!*    4.4    dzs/dy term
!            -----------
!
        imaxj = min(jjss+jnext-1,nsso)
!
        zdzsdy(jiss,jjss:imaxj) = ( xssqo(jiss,inext,il) - xssqo(jiss,jjss,jl)) &
                                    / float(jnext) / zdyeff 
       !
        gdzsdy(jiss,jjss:imaxj) = .true.
!
!
!*    4.5    bottom data point not on the bottom of the grid mesh (one more computation)
!            ----------------------------------------------------
!
        IF (gbound .AND. jjss/=1) THEN
          DO jprev=1,nsso
            IF (jjss-jprev<1) THEN
              il    = ibottom(jl)
              iprev = jjss-jprev+nsso
            ELSE
              il    = jl
              iprev = jjss-jprev
            END IF
            ! no left point
             IF (il==0) EXIT
            ! subgrid data found
            IF (lssqo(jiss,iprev,il)) EXIT
          END DO

          IF (.NOT. (jprev>=nsso+1 .OR. il==0)) THEN
            zdzsdy(jiss,1:jjss) = ( xssqo(jiss,jjss,jl) - xssqo(jiss,iprev,il)) &
                                    / float(jprev) / zdyeff 
!
            gdzsdy(jiss,1:jjss) = .true.
          END IF
        END IF
!
!
        gbound=.false.
!
!
!
!*    4.6   end of loop on jss index
!            ------------------------
!
      END DO
!
!*    4.7    end of loop on iss index
!            ------------------------
!
    END DO
!----------------------------------------------------------------------------
!
!*    5.     Computations of tensor terms
!            ----------------------------
!
!
!*    5.1    test to know if term Hxy is computable
!            --------------------------------------
!
!* 2 values are necessary in the grid point to compute anisotropy
!
    IF ( count(gdzsdx(:,:).AND.gdzsdy(:,:)) ==0 ) cycle
!
!
!*    5.2    SSO quantities are computable
!            -----------------------------

    osso(jl)=.true.
    osso_anis(jl)=count(gdzsdx(:,:).AND.gdzsdy(:,:))>1
!
!
!*    5.3    term Hxx
!            --------
!
    zhxx(jl) = sum(zdzsdx(:,:)*zdzsdx(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))&
                /count(gdzsdx(:,:).AND.gdzsdy(:,:))  
!
!*    5.4    term Hyy
!            --------
!
    zhyy(jl) = sum(zdzsdy(:,:)*zdzsdy(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))&
                /count(gdzsdx(:,:).AND.gdzsdy(:,:)) 
!
!*    5.5    term Hxy
!            --------
!
    zhxy(jl) = sum(zdzsdx(:,:)*zdzsdy(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))&
                /count(gdzsdx(:,:).AND.gdzsdy(:,:))  
!
!-------------------------------------------------------------------------------
!
!*    6.     Next MESONH grid point
!            ----------------------
!
END DO
!
!-------------------------------------------------------------------------------
!
!*    7.     Diagonalization of the tensor
!            -----------------------------
!
WHERE (osso(:))
  zk(:)=0.5*(zhxx(:)+zhyy(:))
  zl(:)=0.5*(zhxx(:)-zhyy(:))
  zm(:)=     zhxy(:)
END WHERE
!
!-------------------------------------------------------------------------------
!
!*    8.     S.S.O. characteristics
!            ----------------------
!
!*    8.1    S.S.O. direction of main axis
!            -----------------------------
!
WHERE (osso(:) .AND. zl(:)>1.e-30 )
  uss%XSSO_DIR(:) = 0.5* atan(zm/zl) * (180./xpi)
END WHERE
!
WHERE (osso(:) .AND. zl(:)<-1.e-30 )
  uss%XSSO_DIR(:) = 0.5* atan(zm/zl) * (180./xpi) + 90.
END WHERE
!
WHERE (osso(:) .AND. abs(zl(:))<=1.e-30 )
  uss%XSSO_DIR(:) = 45.
END WHERE
!
WHERE (osso(:) .AND. uss%XSSO_DIR(:)>90. )
  uss%XSSO_DIR(:) = uss%XSSO_DIR(:) - 180.
END WHERE
!
!*    8.2    S.S.O. slope
!            ------------
!
WHERE (osso(:))
  uss%XSSO_SLOPE(:) = sqrt( zk+sqrt(zl*zl+zm*zm) )
END WHERE
!
!*    8.3    S.S.O. anisotropy
!            -----------------
!
WHERE (osso_anis(:) .AND. (zk+sqrt(zl*zl+zm*zm)) >0. )
  uss%XSSO_ANIS(:)=sqrt( max(zk-sqrt(zl*zl+zm*zm),0.) / (zk+sqrt(zl*zl+zm*zm)))
END WHERE
!
WHERE (osso_anis(:) .AND. (zk+sqrt(zl*zl+zm*zm))==0. )
  uss%XSSO_ANIS(:)=1.
END WHERE
IF (lhook) CALL dr_hook('SSO',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE sso