pgd_topd.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 pgd_topd (HISBA, HGRID, PGRID_PAR, KDIM_FULL, PSSO_SLOPE, HPROGRAM)
!     #############################################################
!
!!****  *PGD_TOPD* - routine to determine the masks that permit to couple ISBA grid with Topmodel one
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!   (from initial version of init_coupl_topo.f90 by K. Chancibault, modified by
!!     M. Lelay and B. Vincendon)
!!
!!    AUTHOR
!!    ------
!!      B. Vincendon   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    11/2011
!!                 03/2014 (E. Artinian) manages the option CGRID='IGN'
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_topd_par, ONLY : nunit
USE modd_topodyn,         ONLY : ccat, nncat, xrtop_d2, nmesht, xdxt
USE modd_coupling_topd,   ONLY : lcoupl_topd, nimax, njmax, &
                                 xxi, xyi, nmaski, nmaskt, nnpix,&
                                 nnbv_in_mesh, xbv_in_mesh, xtotbv_in_mesh
USE modd_dummy_exp_profile, ONLY : xf_param_bv, xc_depth_ratio_bv
!
USE modd_surf_par,          ONLY : nundef
!
!
USE mode_gridtype_conf_proj
USE mode_gridtype_lonlat_reg
USE mode_gridtype_ign
!
USE modi_get_luout
USE modi_read_nam_pgd_topd
USE modi_init_topd_pgd
USE modi_abor1_sfx
USE modi_make_mask_topd_to_isba
USE modi_make_mask_isba_to_topd
USE modi_write_file_masktopd
USE modi_open_file
USE modi_close_file
USE modi_topd_to_isba_slope
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
 CHARACTER(LEN=*), INTENT(IN) :: HISBA
 CHARACTER(LEN=*), INTENT(IN) :: HGRID
REAL, DIMENSION(:), INTENT(IN) :: PGRID_PAR
 INTEGER, INTENT(IN) :: KDIM_FULL
 REAL, DIMENSION(:), INTENT(INOUT) :: PSSO_SLOPE
!
CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM    !
!
CHARACTER(LEN=50),DIMENSION(NNCAT) :: CNAME
INTEGER                   :: IL                     ! number of points
INTEGER                   :: JJ,JI,JK,JWRK ! loop control 
INTEGER                   :: JCAT,JMESH,JPIX ! loop control 
INTEGER                           :: ILUOUT       ! Logical unit for output filr
INTEGER                           :: IMESHL       !  number of ISBA grid nodes
INTEGER                           :: ILAMBERT     ! Lambert projection type
!
REAL, DIMENSION(:), ALLOCATABLE   :: ZXI, ZYI     ! natural coordinates of ISBA grid (conformal projection or latlon)
REAL, DIMENSION(:), ALLOCATABLE   :: ZDXI, ZDYI   ! Isba grid resolution in the conformal projection
REAL, DIMENSION(:), ALLOCATABLE   :: ZXN, ZYN     ! isba nodes coordinates in the Lambert II coordinates
REAL, DIMENSION(:), ALLOCATABLE   :: ZLAT,ZLON    ! Isba nodes geographical coordinates
REAL, DIMENSION(:), ALLOCATABLE   :: ZDTAV        ! Averaged depth soil on TOP-LAT grid
REAL                              :: ZLAT0    ! reference latitude
REAL                              :: ZLON0    ! reference longitude
REAL                              :: ZLONMIN,ZLONMAX  ! min and max longitude values (latlon coordinates)
REAL                              :: ZLATMIN,ZLATMAX  ! min and max latitude values (latlon coordinates)
REAL                              :: ZRPK     ! projection parameter 
!                                             !   K=1 : stereographic north pole
!                                             ! 0<K<1 : Lambert, north hemisphere
!                                             !   K=0 : Mercator
!                                             !-1<K<0 : Lambert, south hemisphere
!                                             !   K=-1: stereographic south pole
REAL                              :: ZBETA    ! angle between grid and reference longitude
REAL                              :: ZLATOR   ! latitude  of point of coordinates X=0, Y=0
REAL                              :: ZLONOR   ! longitude of point of coordinates X=0, Y=0!
REAL, DIMENSION(:), ALLOCATABLE   :: ZF_PARAM,ZC_DEPTH_RATIO  !
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('PGD_TOPD',0,zhook_handle)
!
CALL GET_LUOUT(HPROGRAM,ILUOUT)
!  
CALL READ_NAM_PGD_TOPD(HPROGRAM,LCOUPL_TOPD,CCAT,XF_PARAM_BV,XC_DEPTH_RATIO_BV)
!
IF (lcoupl_topd .AND. (hisba/='3-L'.AND. hisba/='DIF')) &
  CALL abor1_sfx("PGD_TOPD: coupling with topmodel only runs with CISBA=3-L or CISBA=DIF  ")
!
!         1.   Reads the namelists
!              --------------------
IF (lcoupl_topd) THEN
  !
  WRITE(iluout,*) 'Debut pgd_topd'
  !
  !         3.   Initialises variables specific to Topmodel
  !              -------------------------------------------
  WRITE(iluout,*) 'NNCAT',nncat
  !
  CALL init_topd_pgd(hprogram)
  !
  !         4.   Compute masks to couple ISBA and TOPODYN grids
  !              -------------------------------------------------------
  !
  !*      1.      Calcul of the plane coordinates 
  !               -------------------------------
  !
  !*      1.1 Gestion des projections conformes
  !           ---------------------------------
  !
  ALLOCATE(nmaskt(nncat,nmesht))
  nmaskt(:,:)=nundef
  !
  IF(hgrid.EQ.'CONF PROJ') THEN
    !
    WRITE(iluout,*) 'GRILLE PROJ CONF (application Cevennes)'
    !
    !*      1.1.1   lecture des coordonnees X et Y conformes 
    !               -----------------------------------------
    !
    ALLOCATE(zxi(kdim_full))
    ALLOCATE(zyi(kdim_full))
    zxi(:)=0.0
    zyi(:)=0.0
    !
    ALLOCATE(zdxi(kdim_full))
    ALLOCATE(zdyi(kdim_full))
    !
    CALL get_gridtype_conf_proj(pgrid_par,plat0=zlat0,plon0=zlon0,prpk=zrpk, &
                                pbeta=zbeta,plator=zlator,plonor=zlonor,     &
                                kimax=nimax,kjmax=njmax,px=zxi,py=zyi,       &
                                pdx=zdxi,pdy=zdyi)
    !
    imeshl = (nimax+1)*(njmax+1)
    !
    !*      1.1.2   calcul des coordonnees X et Y conformes des noeuds de la grille ISBA
    !               --------------------------------------------------------------------
    ALLOCATE(zxn(imeshl))
    ALLOCATE(zyn(imeshl))
    !
    DO jj=1,nimax
      zxn(jj) = zxi(jj) - zdxi(jj)/2.
    ENDDO
    zxn(nimax+1) = zxi(nimax) + zdxi(nimax)/2.
    !
    DO jj=1,njmax
      jwrk = (jj-1)*(nimax+1)+1   ! indice sur la grille des noeuds
      ji = (jj-1)*nimax+1         ! indice sur la grille des mailles
      zyn(jwrk) = zyi(ji) - zdyi(ji)/2.
    ENDDO
    !
    jj = ((njmax+1)-1)*(nimax+1)+1  ! Indice sur la grille des noeuds
    ji = (njmax-1)*nimax+1
    zyn(jj) = zyi(ji) + zdyi(ji)/2.
    !
    DEALLOCATE(zdxi)
    DEALLOCATE(zdyi)
    DEALLOCATE(zxi)
    DEALLOCATE(zyi)
    !
    DO jj=1,nimax+1
      DO ji=2,njmax+1
        jk = (ji-1)*(nimax+1)+jj
        zxn(jk) = zxn(jj)
      ENDDO
    ENDDO
    !
    DO ji=1,njmax+1
      DO jj=2,nimax+1
        jk = (ji-1)*(nimax+1)+jj
        jwrk = (ji-1)*(nimax+1)+1
        zyn(jk) = zyn(jwrk)
      ENDDO
    ENDDO
    !    
    !*      1.1.3   calcul des coordonnées géographiques des noeuds de la grille ISBA
    !               -----------------------------------------------------------------
    ALLOCATE(zlat(imeshl))
    ALLOCATE(zlon(imeshl))
    CALL latlon_conf_proj(zlat0,zlon0,zrpk,zbeta,zlator,zlonor,zxn,zyn,zlat,zlon)
    DEALLOCATE(zxn)
    DEALLOCATE(zyn)
    !
    !*      1.2 Gestion des coordonnees geographiques
    !           -------------------------------------
    !
  ELSE IF(hgrid.EQ.'LONLAT REG') THEN
    !
    WRITE(iluout,*) 'GRILLE LONLAT REG (application AMMA)' 
    !
    ALLOCATE(zxi(kdim_full))
    ALLOCATE(zyi(kdim_full))
    zxi(:)=0.0
    zyi(:)=0.0
    CALL get_gridtype_lonlat_reg(pgrid_par,plonmin=zlonmin,plonmax=zlonmax,             &
                                 platmin=zlatmin,platmax=zlatmax,klon=nimax,klat=njmax, &
                                 kl=il,plon=zxi,plat=zyi)
    !
    imeshl=(nimax+1)*(njmax+1)
    !
    ALLOCATE(zlon(imeshl))
    ALLOCATE(zlat(imeshl))
    ALLOCATE(zdxi(kdim_full))
    ALLOCATE(zdyi(kdim_full))
    !
    zdxi(:)=(zlonmax-zlonmin)/(nimax-1)
    zdyi(:)=(zlatmax-zlatmin)/(njmax-1)
    !
    DO jj=1,nimax
      zlon(jj) = zxi(jj) - zdxi(jj)/2.
    ENDDO
    zlon(nimax+1) = zxi(nimax) + zdxi(nimax)/2.
    !
    DO jj=1,njmax
      jwrk=(jj-1)*(nimax+1)+1   ! indice sur la grille des noeuds
      ji=(jj-1)*nimax+1        ! indice sur la grille des mailles
      zlat(jwrk) = zyi(ji) - zdyi(ji)/2.
    ENDDO
    !
    jj=((njmax+1)-1)*(nimax+1)+1  ! Indice sur la grille des noeuds
    ji=(njmax-1)*nimax+1
    zlat(jj) = zyi(ji) + zdyi(ji)/2.
    !
    DEALLOCATE(zdxi)
    DEALLOCATE(zdyi)
    DEALLOCATE(zxi)
    DEALLOCATE(zyi)
    !
    DO jj=1,nimax+1
      DO ji=2,njmax+1
        jk = (ji-1)*(nimax+1)+jj
        zlon(jk) = zlon(jj)
      ENDDO
    ENDDO
    !
    DO ji=1,njmax+1
      DO jj=2,nimax+1
        jk=(ji-1)*(nimax+1)+jj
        jwrk=(ji-1)*(nimax+1)+1
        zlat(jk)=zlat(jwrk)
      ENDDO
    ENDDO
    !
  ! Modification by Eram Artinian to take into account IGN grid 1
  ELSE IF (hgrid=='IGN') THEN 
    WRITE(iluout,*) 'GRILLE IGN (application Bulgarie)' 
    ALLOCATE(zxn(kdim_full))
    ALLOCATE(zyn(kdim_full))
    CALL get_gridtype_ign(pgrid_par,klambert=ilambert,&
                          kl=il,px=zxn,py=zyn,kdimx=nimax)
    imeshl=il
    ALLOCATE(zlat(imeshl))
    ALLOCATE(zlon(imeshl))
    CALL latlon_ign(ilambert,zxn,zyn,zlat,zlon)
  !End modification by Eram Artinian to take into account IGN grid
  ELSE 
    !       
    WRITE(iluout,*) 'ERREUR: TYPE DE GRILLE NON GERE PAR LE CODE'
    CALL abor1_sfx("PGD_TOPD: TYPE DE GRILLE NON GERE PAR LE CODE")
    !
  ENDIF
  !
  !*      2.0  calcul des coordonnées lambert II étendu des noeuds de la grille ISBA
  !            ----------------------------------------------------------------------
  !
  ALLOCATE(xxi(imeshl))
  ALLOCATE(xyi(imeshl))
  ! Modification by Eram Artinian to take into account IGN grid 2
  IF (hgrid/='IGN') THEN
    CALL xy_ign(5,xxi,xyi,zlat,zlon)
  ELSE
    xxi=zxn
    xyi=zyn
    DEALLOCATE(zxn)
    DEALLOCATE(zyn)
  ENDIF
  !End modification by Eram Artinian to take into account IGN grid 2
  DEALLOCATE(zlat)
  DEALLOCATE(zlon)
  !
  !*      2.0     mask
  !               ----
  !***
  CALL make_mask_topd_to_isba(hgrid, pgrid_par, kdim_full)
  !
  ALLOCATE(nnpix(kdim_full))
  nnpix(:) = nundef
  DO jj=1,kdim_full
    nnpix(jj) = count(nmaskt(:,:)==jj)
  ENDDO
  !
  CALL make_mask_isba_to_topd(kdim_full)
  !
  CALL write_file_masktopd(kdim_full)
  !
  !*        3.0 Compute Mean slope over each ISBA_MESH
!            ----------------------------------------------------------------------
   CALL topd_to_isba_slope(psso_slope, kdim_full)
!
!*        4.0  Compute F and DC for each ISBA mesh
  !            ----------------------------------------------------------------------
  !
  ALLOCATE(nnbv_in_mesh(kdim_full,nncat))
  ALLOCATE(xbv_in_mesh(kdim_full,nncat))
  ALLOCATE(xtotbv_in_mesh(kdim_full))
  !
  xtotbv_in_mesh(:) = 0.0
  !
  DO jmesh=1,kdim_full
    xbv_in_mesh(jmesh,:)=0.0
    DO jcat=1,nncat
      nnbv_in_mesh(jmesh,jcat) = count(nmaski(jmesh,jcat,:)/=nundef)
      xbv_in_mesh(jmesh,jcat) = REAL(nnbv_in_mesh(jmesh,jcat))*XDXT(jcat)**2
      xtotbv_in_mesh(jmesh) = xtotbv_in_mesh(jmesh) + xbv_in_mesh(jmesh,jcat)
    ENDDO
  ENDDO
  !
  ALLOCATE (zf_param(kdim_full))
  ALLOCATE (zc_depth_ratio(kdim_full))
  !
  zf_param(:) = 0.
  zc_depth_ratio(:) = 0.
  DO jcat=1,nncat
    DO jmesh=1,kdim_full
      IF ( xtotbv_in_mesh(jmesh)/=0. ) THEN
        zf_param(jmesh) = zf_param(jmesh) + xf_param_bv(jcat)*xbv_in_mesh(jmesh,jcat)/xtotbv_in_mesh(jmesh)
        zc_depth_ratio(jmesh) = zc_depth_ratio(jmesh) + xc_depth_ratio_bv(jcat)*xbv_in_mesh(jmesh,jcat)/xtotbv_in_mesh(jmesh)
      ENDIF
    ENDDO
  ENDDO
  !
  WHERE (zf_param==0.)
    zf_param=2.5
    zc_depth_ratio=1.
  ENDWHERE
  !
  !write(*,*) 'f min max isba',MINVAL(ZF_PARAM),MAXVAL(ZF_PARAM)
  !write(*,*) 'dc min max isba',MINVAL(ZC_DEPTH_RATIO),MAXVAL(ZC_DEPTH_RATIO)
  !
  CALL open_file('ASCII ',nunit,'carte_f_dc.txt','FORMATTED',haction='WRITE')
  DO jmesh=1,kdim_full
    WRITE(nunit,*) zf_param(jmesh),zc_depth_ratio(jmesh)
  ENDDO
  CALL close_file('ASCII ',nunit)
  !
  DEALLOCATE(zf_param)
  DEALLOCATE(zc_depth_ratio)
  !
  WRITE(iluout,*) 'Couplage avec TOPMODEL active'
  !
ELSE
  !
  WRITE(iluout,*) 'Pas de couplage avec TOPMODEL'
  !
ENDIF
!
IF (lhook) CALL dr_hook('PGD_TOPD',1,zhook_handle)
!
END SUBROUTINE pgd_topd