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 (U, UG, USS, OSSO, OSSO_ANIS)
!     #########################
!
!!*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_surfex_mpi, ONLY : nrank
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_sso_n, ONLY : sso_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_surfex_mpi, ONLY : nrank, npio
!
USE modi_get_mesh_dim
USE modi_get_adjacent_meshes
!
USE modi_gather_and_write_mpi
USE modi_read_and_send_mpi
!
USE modd_csts,           ONLY : xpi
USE modd_pgdwork,        ONLY : nsso, xssqo, lssqo
USE modd_pgd_grid,       ONLY : nl, cgrid
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of dummy arguments
!            ------------------------------
!
!
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(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
!
!*    0.2    Declaration of indexes
!            ----------------------
!
INTEGER :: JS1, JS2
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(:), ALLOCATABLE :: ILEFT, IRIGHT, ITOP, IBOTTOM
INTEGER :: ICOUNT
!
!*    0.3    Declaration of working arrays inside a MESONH grid (JI,JJ)
!            -----------------------------
!
INTEGER, DIMENSION(:), ALLOCATABLE :: ISSO, ISSOT
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: ISSQOT, ISSQO
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSSQO
LOGICAL, DIMENSION(:,:,:), ALLOCATABLE :: GSSQO
LOGICAL, DIMENSION(:), ALLOCATABLE :: GSSO, GSSO_ANIS
!
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(:), ALLOCATABLE :: ZDX, ZDY, ZSEA, ZMESH_SIZE
REAL, DIMENSION(:), ALLOCATABLE :: ZHXX0, ZHXY0, ZHYY0
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)
!
!*    1.1    Occurence of computation of the coefficients
!            --------------------------------------------
!
osso(:)=.false.
osso_anis(:)=.false.
!
zhyy(:) = 0.
zhxx(:) = 0.
zhxy(:) = 0.
!
!*    1.2    Grid dimension (meters)
!            -----------------------
!
IF (nrank==npio) THEN
  ALLOCATE(zsea(u%NDIM_FULL))
  ALLOCATE(zmesh_size(u%NDIM_FULL))
  ALLOCATE(zssqo(u%NDIM_FULL,nsso,nsso))
  ALLOCATE(issqot(u%NDIM_FULL,nsso,nsso))
ELSE
  ALLOCATE(zsea(0))
  ALLOCATE(zssqo(0,0,0))
  ALLOCATE(issqot(0,0,0))
ENDIF
!
 CALL gather_and_write_mpi(u%XSEA,zsea)
 CALL gather_and_write_mpi(ug%G%XMESH_SIZE,zmesh_size) 
 CALL gather_and_write_mpi(xssqo,zssqo)
!
ALLOCATE(issqo(SIZE(xssqo,1),nsso,nsso))
IF (SIZE(issqo)/=0) THEN
  issqo(:,:,:) = 0
  WHERE (lssqo(:,:,:)) issqo(:,:,:) = 1
ENDIF
 CALL gather_and_write_mpi(issqo,issqot)
DEALLOCATE(issqo)
!
IF (nrank==npio) THEN
  !
  ALLOCATE(zdx(u%NDIM_FULL),zdy(u%NDIM_FULL))
  CALL get_mesh_dim(cgrid,ug%NGRID_FULL_PAR,u%NDIM_FULL,ug%XGRID_FULL_PAR,zdx,zdy,zmesh_size)
  DEALLOCATE(zmesh_size)

  !
  !
  !*    1.3    Left, top, right and bottom adjacent gris meshes
  !            ------------------------------------------------
  !
  ALLOCATE(ileft(u%NDIM_FULL),iright(u%NDIM_FULL),itop(u%NDIM_FULL),ibottom(u%NDIM_FULL))
  CALL get_adjacent_meshes(cgrid,ug%NGRID_FULL_PAR,u%NDIM_FULL,ug%XGRID_FULL_PAR,&
                           ileft,iright,itop,ibottom)
  !
  ALLOCATE(gsso(u%NDIM_FULL),gsso_anis(u%NDIM_FULL))
  ALLOCATE(zhxx0(u%NDIM_FULL),zhxy0(u%NDIM_FULL),zhyy0(u%NDIM_FULL))
  gsso(:) = .false.
  gsso_anis(:) = .false.
  zhyy0(:) = 0.
  zhxx0(:) = 0.
  zhxy0(:) = 0.
  !
  !
  !*    2.     Loop on MESONH grid points
  !            --------------------------
  !
  !
  !done on the whole grid (needed because of use of near points)
  DO jl=1,u%NDIM_FULL
    !
    !
    !*    2.1    No land in grid mesh
    !            --------------------
    !
    IF (zsea(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( issqot(jl,:,:)==1) == 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 ( issqot(jl,jiss,jjss)==0 ) 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 (issqot(il,inext,jjss)==1) 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) = ( zssqo(il,inext,jjss) - zssqo(jl,jiss,jjss)) &
                                    / 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 (issqot(il,iprev,jjss)==1) EXIT
          END DO

          IF (.NOT. (jprev>=nsso+1 .OR. il==0)) THEN
            zdzsdx(1:jiss,jjss) = ( zssqo(jl,jiss,jjss) - zssqo(il,iprev,jjss)) &
                                    / 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 (issqot(jl,jiss,jjss)==0 ) 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 (issqot(il,jiss,inext)==1) 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) = ( zssqo(il,jiss,inext) - zssqo(jl,jiss,jjss)) &
                                    / 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 (issqot(il,jiss,iprev)==1) EXIT
          END DO

          IF (.NOT. (jprev>=nsso+1 .OR. il==0)) THEN
            zdzsdy(jiss,1:jjss) = ( zssqo(jl,jiss,jjss) - zssqo(il,jiss,iprev)) &
                                    / 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
    !            --------------------------------------
    !
    icount =  count(gdzsdx(:,:).AND.gdzsdy(:,:))
    !* 2 values are necessary in the grid point to compute anisotropy
    !
    IF ( icount==0 ) cycle
    !
    !
    !*    5.2    SSO quantities are computable
    !            -----------------------------
    !
    gsso(jl)=.true.
    gsso_anis(jl)=icount>1
    !
    !
    !*    5.3    term Hxx
    !            --------
    !
    zhxx0(jl) = sum(zdzsdx(:,:)*zdzsdx(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))/icount  
    !
    !*    5.4    term Hyy
    !            --------
    !
    zhyy0(jl) = sum(zdzsdy(:,:)*zdzsdy(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))/icount 
    !
    !*    5.5    term Hxy
    !            --------
    !
    zhxy0(jl) = sum(zdzsdx(:,:)*zdzsdy(:,:),mask=gdzsdx(:,:).AND.gdzsdy(:,:))/icount  
    !
    !-------------------------------------------------------------------------------
    !
    !*    6.     Next MESONH grid point
    !            ----------------------
    !
  END DO    
  !
ELSE
  !
  ALLOCATE(zhxx0(0),zhyy0(0),zhxy0(0))
  !
ENDIF    
!
DEALLOCATE(zssqo,issqot)
DEALLOCATE(zsea)
!
!each contrib sent to each task
 CALL read_and_send_mpi(zhxx0,zhxx)
 CALL read_and_send_mpi(zhyy0,zhyy)
 CALL read_and_send_mpi(zhxy0,zhxy)
DEALLOCATE(zhxx0,zhyy0,zhxy0)
!
IF (nrank==npio) THEN
  ALLOCATE(issot(u%NDIM_FULL))
  issot(:) = 0
  WHERE (gsso(:)) issot(:) = 1
  DEALLOCATE(gsso)
ELSE
  ALLOCATE(issot(0))
ENDIF
ALLOCATE(isso(nl))
 CALL read_and_send_mpi(issot,isso)
WHERE(isso(:)==1) osso(:) = .true.
!
IF (nrank==npio) THEN
  issot(:) = 0
  WHERE (gsso_anis(:)) issot(:) = 1
  DEALLOCATE(gsso_anis)
ENDIF
 CALL read_and_send_mpi(issot,isso)
WHERE(isso(:)==1) osso_anis(:) = .true.
DEALLOCATE(isso)
!
DEALLOCATE(issot)
!
!----------------------------------------------------------------------------
!
!*    7.     Diagonalization of the tensor
!     -----------------------------
!
zk=0.
zl=0.
zm=0.
!
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