SURFEX v7.3
General documentation of Surfex
 All Classes Files Functions Variables Typedefs
/home/dasprezs/EXPORT_v7_3/src/LIB/TOPD/pgd_topd.F90
Go to the documentation of this file.
00001 !-------------------------------------------------------------------------------
00002 !     #############################################################
00003       SUBROUTINE PGD_TOPD(HPROGRAM)
00004 !     #############################################################
00005 !
00006 !!****  *PGD_TOPD* - routine to determine the masks that permit to couple ISBA grid with Topmodel one
00007 !!
00008 !!    PURPOSE
00009 !!    -------
00010 !!
00011 !!**  METHOD
00012 !!    ------
00013 !!
00014 !!    EXTERNAL
00015 !!    --------
00016 !!
00017 !!    IMPLICIT ARGUMENTS
00018 !!    ------------------
00019 !!
00020 !!    REFERENCE
00021 !!    ---------
00022 !!   (from initial version of init_coupl_topo.f90 by K. Chancibault, modified by
00023 !!     M. Lelay and B. Vincendon)
00024 !!
00025 !!    AUTHOR
00026 !!    ------
00027 !!      B. Vincendon   *Meteo France*   
00028 !!
00029 !!    MODIFICATIONS
00030 !!    -------------
00031 !!      Original    11/2011
00032 !-------------------------------------------------------------------------------
00033 !
00034 !*       0.    DECLARATIONS
00035 !              ------------
00036 !
00037 USE MODD_TOPODYN,         ONLY : CCAT, NNCAT, XRTOP_D2, NMESHT, XDXT
00038 USE MODD_COUPLING_TOPD,   ONLY : LCOUPL_TOPD, NIMAX, NJMAX, &
00039                                  XXI, XYI, NMASKI, NMASKT, NNPIX,&
00040                                  NNBV_IN_MESH, XBV_IN_MESH, XTOTBV_IN_MESH
00041 USE MODD_DUMMY_EXP_PROFILE, ONLY : XF_PARAM_BV, XC_DEPTH_RATIO_BV
00042 !
00043 USE MODD_SURF_PAR,          ONLY : NUNDEF
00044 USE MODD_SURF_ATM_GRID_N,   ONLY : XGRID_PAR, CGRID  !mll: ajout CGRID
00045 USE MODD_SURF_ATM_n,        ONLY : NDIM_FULL
00046 !
00047 USE MODD_ISBA_n, ONLY : CISBA
00048 !
00049 USE MODE_GRIDTYPE_CONF_PROJ
00050 USE MODE_GRIDTYPE_LONLAT_REG
00051 USE MODE_GRIDTYPE_IGN
00052 !
00053 USE MODI_GET_LUOUT
00054 USE MODI_READ_NAM_PGD_TOPD
00055 USE MODI_INIT_TOPD
00056 USE MODI_ABOR1_SFX
00057 USE MODI_MAKE_MASK_TOPD_TO_ISBA
00058 USE MODI_MAKE_MASK_ISBA_TO_TOPD
00059 USE MODI_WRITE_FILE_MASKTOPD
00060 USE MODI_OPEN_FILE
00061 USE MODI_CLOSE_FILE
00062 USE MODI_TOPD_TO_ISBA_SLOPE
00063 !
00064 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00065 USE PARKIND1  ,ONLY : JPRB
00066 !
00067 IMPLICIT NONE
00068 !
00069 !*       0.1   Declarations of arguments
00070 !              -------------------------
00071 !
00072  CHARACTER(LEN=*),  INTENT(IN)     :: HPROGRAM    !
00073 !
00074  CHARACTER(LEN=50),DIMENSION(NNCAT) :: CNAME
00075 INTEGER                   :: IL                     ! number of points
00076 INTEGER                   :: JJ,JI,JK,JWRK ! loop control 
00077 INTEGER                   :: JCAT,JMESH,JPIX ! loop control 
00078 INTEGER                           :: ILUOUT       ! Logical unit for output filr
00079 INTEGER                           :: IUNIT       ! Logical unit for carte_asat file
00080 INTEGER                           :: IMESHL       !  number of ISBA grid nodes
00081 !
00082 REAL, DIMENSION(:), ALLOCATABLE   :: ZXI, ZYI     ! natural coordinates of ISBA grid (conformal projection or latlon)
00083 REAL, DIMENSION(:), ALLOCATABLE   :: ZDXI, ZDYI   ! Isba grid resolution in the conformal projection
00084 REAL, DIMENSION(:), ALLOCATABLE   :: ZXN, ZYN     ! isba nodes coordinates in the Lambert II coordinates
00085 REAL, DIMENSION(:), ALLOCATABLE   :: ZLAT,ZLON    ! Isba nodes geographical coordinates
00086 REAL, DIMENSION(:), ALLOCATABLE   :: ZDTAV        ! Averaged depth soil on TOP-LAT grid
00087 REAL                              :: ZLAT0    ! reference latitude
00088 REAL                              :: ZLON0    ! reference longitude
00089 REAL                              :: ZLONMIN,ZLONMAX  ! min and max longitude values (latlon coordinates)
00090 REAL                              :: ZLATMIN,ZLATMAX  ! min and max latitude values (latlon coordinates)
00091 REAL                              :: ZRPK     ! projection parameter 
00092 !                                             !   K=1 : stereographic north pole
00093 !                                             ! 0<K<1 : Lambert, north hemisphere
00094 !                                             !   K=0 : Mercator
00095 !                                             !-1<K<0 : Lambert, south hemisphere
00096 !                                             !   K=-1: stereographic south pole
00097 REAL                              :: ZBETA    ! angle between grid and reference longitude
00098 REAL                              :: ZLATOR   ! latitude  of point of coordinates X=0, Y=0
00099 REAL                              :: ZLONOR   ! longitude of point of coordinates X=0, Y=0!
00100 REAL, DIMENSION(:), ALLOCATABLE   :: ZF_PARAM,ZC_DEPTH_RATIO  !
00101 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00102 !-------------------------------------------------------------------------------
00103 IF (LHOOK) CALL DR_HOOK('PGD_TOPD',0,ZHOOK_HANDLE)
00104 !
00105  CALL GET_LUOUT(HPROGRAM,ILUOUT)
00106 !  
00107  CALL READ_NAM_PGD_TOPD(HPROGRAM,LCOUPL_TOPD,CCAT,XF_PARAM_BV,XC_DEPTH_RATIO_BV)
00108 !
00109 IF (LCOUPL_TOPD .AND. CISBA/='3-L') &
00110   CALL ABOR1_SFX("PGD_TOPD: coupling with topmodel only runs with CISBA=3-L")
00111 !
00112 !         1.   Reads the namelists
00113 !              --------------------
00114 IF (LCOUPL_TOPD) THEN
00115   !
00116   WRITE(ILUOUT,*) 'Debut pgd_topd'
00117   !
00118   !         3.   Initialises variables specific to Topmodel
00119   !              -------------------------------------------
00120   WRITE(ILUOUT,*) 'NNCAT',NNCAT
00121   !
00122   CALL INIT_TOPD(HPROGRAM)
00123   !
00124   !         4.   Compute masks to couple ISBA and TOPODYN grids
00125   !              -------------------------------------------------------
00126   !
00127   !*      1.      Calcul of the plane coordinates 
00128   !               -------------------------------
00129   !
00130   !*      1.1 Gestion des projections conformes
00131   !           ---------------------------------
00132   !
00133   ALLOCATE(NMASKT(NNCAT,NMESHT))
00134   NMASKT(:,:)=NUNDEF
00135   !
00136   IF(CGRID.EQ.'CONF PROJ') THEN
00137     !
00138     WRITE(ILUOUT,*) 'GRILLE PROJ CONF (application Cevennes)'
00139     !
00140     !*      1.1.1   lecture des coordonnees X et Y conformes 
00141     !               -----------------------------------------
00142     !
00143     ALLOCATE(ZXI(NDIM_FULL))
00144     ALLOCATE(ZYI(NDIM_FULL))
00145     ZXI(:)=0.0
00146     ZYI(:)=0.0
00147     !
00148     ALLOCATE(ZDXI(NDIM_FULL))
00149     ALLOCATE(ZDYI(NDIM_FULL))
00150     !
00151     CALL GET_GRIDTYPE_CONF_PROJ(XGRID_PAR,PLAT0=ZLAT0,PLON0=ZLON0,PRPK=ZRPK, &
00152                                 PBETA=ZBETA,PLATOR=ZLATOR,PLONOR=ZLONOR,     &
00153                                 KIMAX=NIMAX,KJMAX=NJMAX,PX=ZXI,PY=ZYI,       &
00154                                 PDX=ZDXI,PDY=ZDYI)
00155     !
00156     IMESHL = (NIMAX+1)*(NJMAX+1)
00157     !
00158     !*      1.1.2   calcul des coordonnees X et Y conformes des noeuds de la grille ISBA
00159     !               --------------------------------------------------------------------
00160     ALLOCATE(ZXN(IMESHL))
00161     ALLOCATE(ZYN(IMESHL))
00162     !
00163     DO JJ=1,NIMAX
00164       ZXN(JJ) = ZXI(JJ) - ZDXI(JJ)/2.
00165     ENDDO
00166     ZXN(NIMAX+1) = ZXI(NIMAX) + ZDXI(NIMAX)/2.
00167     !
00168     DO JJ=1,NJMAX
00169       JWRK = (JJ-1)*(NIMAX+1)+1   ! indice sur la grille des noeuds
00170       JI = (JJ-1)*NIMAX+1         ! indice sur la grille des mailles
00171       ZYN(JWRK) = ZYI(JI) - ZDYI(JI)/2.
00172     ENDDO
00173     !
00174     JJ = ((NJMAX+1)-1)*(NIMAX+1)+1  ! Indice sur la grille des noeuds
00175     JI = (NJMAX-1)*NIMAX+1
00176     ZYN(JJ) = ZYI(JI) + ZDYI(JI)/2.
00177     !
00178     DEALLOCATE(ZDXI)
00179     DEALLOCATE(ZDYI)
00180     DEALLOCATE(ZXI)
00181     DEALLOCATE(ZYI)
00182     !
00183     DO JJ=1,NIMAX+1
00184       DO JI=2,NJMAX+1
00185         JK = (JI-1)*(NIMAX+1)+JJ
00186         ZXN(JK) = ZXN(JJ)
00187       ENDDO
00188     ENDDO
00189     !
00190     DO JI=1,NJMAX+1
00191       DO JJ=2,NIMAX+1
00192         JK = (JI-1)*(NIMAX+1)+JJ
00193         JWRK = (JI-1)*(NIMAX+1)+1
00194         ZYN(JK) = ZYN(JWRK)
00195       ENDDO
00196     ENDDO
00197     !    
00198     !*      1.1.3   calcul des coordonnées géographiques des noeuds de la grille ISBA
00199     !               -----------------------------------------------------------------
00200     ALLOCATE(ZLAT(IMESHL))
00201     ALLOCATE(ZLON(IMESHL))
00202     CALL LATLON_CONF_PROJ(ZLAT0,ZLON0,ZRPK,ZBETA,ZLATOR,ZLONOR,ZXN,ZYN,ZLAT,ZLON)
00203     DEALLOCATE(ZXN)
00204     DEALLOCATE(ZYN)
00205     !
00206     !*      1.2 Gestion des coordonnees geographiques
00207     !           -------------------------------------
00208     !
00209   ELSE IF(CGRID.EQ.'LONLAT REG') THEN
00210     !
00211     WRITE(ILUOUT,*) 'GRILLE LONLAT REG (application AMMA)' 
00212     !
00213     ALLOCATE(ZXI(NDIM_FULL))
00214     ALLOCATE(ZYI(NDIM_FULL))
00215     ZXI(:)=0.0
00216     ZYI(:)=0.0
00217     CALL GET_GRIDTYPE_LONLAT_REG(XGRID_PAR,PLONMIN=ZLONMIN,PLONMAX=ZLONMAX,             &
00218                                  PLATMIN=ZLATMIN,PLATMAX=ZLATMAX,KLON=NIMAX,KLAT=NJMAX, &
00219                                  KL=IL,PLON=ZXI,PLAT=ZYI)
00220     !
00221     IMESHL=(NIMAX+1)*(NJMAX+1)
00222     !
00223     ALLOCATE(ZLON(IMESHL))
00224     ALLOCATE(ZLAT(IMESHL))
00225     ALLOCATE(ZDXI(NDIM_FULL))
00226     ALLOCATE(ZDYI(NDIM_FULL))
00227     !
00228     ZDXI(:)=(ZLONMAX-ZLONMIN)/(NIMAX-1)
00229     ZDYI(:)=(ZLATMAX-ZLATMIN)/(NJMAX-1)
00230     !
00231     DO JJ=1,NIMAX
00232       ZLON(JJ) = ZXI(JJ) - ZDXI(JJ)/2.
00233     ENDDO
00234     ZLON(NIMAX+1) = ZXI(NIMAX) + ZDXI(NIMAX)/2.
00235     !
00236     DO JJ=1,NJMAX
00237       JWRK=(JJ-1)*(NIMAX+1)+1   ! indice sur la grille des noeuds
00238       JI=(JJ-1)*NIMAX+1        ! indice sur la grille des mailles
00239       ZLAT(JWRK) = ZYI(JI) - ZDYI(JI)/2.
00240     ENDDO
00241     !
00242     JJ=((NJMAX+1)-1)*(NIMAX+1)+1  ! Indice sur la grille des noeuds
00243     JI=(NJMAX-1)*NIMAX+1
00244     ZLAT(JJ) = ZYI(JI) + ZDYI(JI)/2.
00245     !
00246     DEALLOCATE(ZDXI)
00247     DEALLOCATE(ZDYI)
00248     DEALLOCATE(ZXI)
00249     DEALLOCATE(ZYI)
00250     !
00251     DO JJ=1,NIMAX+1
00252       DO JI=2,NJMAX+1
00253         JK = (JI-1)*(NIMAX+1)+JJ
00254         ZLON(JK) = ZLON(JJ)
00255       ENDDO
00256     ENDDO
00257     !
00258     DO JI=1,NJMAX+1
00259       DO JJ=2,NIMAX+1
00260         JK=(JI-1)*(NIMAX+1)+JJ
00261         JWRK=(JI-1)*(NIMAX+1)+1
00262         ZLAT(JK)=ZLAT(JWRK)
00263       ENDDO
00264     ENDDO
00265     !
00266   ELSE
00267     !       
00268     WRITE(ILUOUT,*) 'ERREUR: TYPE DE GRILLE NON GERE PAR LE CODE'
00269     CALL ABOR1_SFX("PGD_TOPD: TYPE DE GRILLE NON GERE PAR LE CODE")
00270     !
00271   ENDIF
00272   !
00273   !*      2.0  calcul des coordonnées lambert II étendu des noeuds de la grille ISBA
00274   !            ----------------------------------------------------------------------
00275   !
00276   ALLOCATE(XXI(IMESHL))
00277   ALLOCATE(XYI(IMESHL))
00278   CALL XY_IGN(5,XXI,XYI,ZLAT,ZLON)
00279   DEALLOCATE(ZLAT)
00280   DEALLOCATE(ZLON)
00281   !
00282   !*      2.0     mask
00283   !               ----
00284   !***
00285   CALL MAKE_MASK_TOPD_TO_ISBA(NDIM_FULL)
00286   !
00287   ALLOCATE(NNPIX(NDIM_FULL))
00288   NNPIX(:) = NUNDEF
00289   DO JJ=1,NDIM_FULL
00290     NNPIX(JJ) = COUNT(NMASKT(:,:)==JJ)
00291   ENDDO
00292   !
00293   CALL MAKE_MASK_ISBA_TO_TOPD(NDIM_FULL)
00294   !
00295   CALL WRITE_FILE_MASKTOPD(NDIM_FULL)
00296   !
00297   !*        3.0 Compute Mean slope over each ISBA_MESH
00298   !            ----------------------------------------------------------------------
00299   CALL TOPD_TO_ISBA_SLOPE(NDIM_FULL)
00300   !
00301   !*        4.0  Compute F and DC for each ISBA mesh
00302   !            ----------------------------------------------------------------------
00303   !
00304   ALLOCATE(NNBV_IN_MESH(NDIM_FULL,NNCAT))
00305   ALLOCATE(XBV_IN_MESH(NDIM_FULL,NNCAT))
00306   ALLOCATE(XTOTBV_IN_MESH(NDIM_FULL))
00307   !
00308   XTOTBV_IN_MESH(:) = 0.0
00309   !
00310   DO JMESH=1,NDIM_FULL
00311     XBV_IN_MESH(JMESH,:)=0.0
00312     DO JCAT=1,NNCAT
00313       NNBV_IN_MESH(JMESH,JCAT) = COUNT(NMASKI(JMESH,JCAT,:)/=NUNDEF)
00314       XBV_IN_MESH(JMESH,JCAT) = REAL(NNBV_IN_MESH(JMESH,JCAT))*XDXT(JCAT)**2
00315       XTOTBV_IN_MESH(JMESH) = XTOTBV_IN_MESH(JMESH) + XBV_IN_MESH(JMESH,JCAT)
00316     ENDDO
00317   ENDDO
00318   !
00319   ALLOCATE (ZF_PARAM(NDIM_FULL))
00320   ALLOCATE (ZC_DEPTH_RATIO(NDIM_FULL))
00321   !
00322   ZF_PARAM(:) = 0.
00323   ZC_DEPTH_RATIO(:) = 0.
00324   DO JCAT=1,NNCAT
00325     DO JMESH=1,NDIM_FULL
00326       IF ( XTOTBV_IN_MESH(JMESH)/=0. ) THEN
00327         ZF_PARAM(JMESH) = ZF_PARAM(JMESH) + XF_PARAM_BV(JCAT)*XBV_IN_MESH(JMESH,JCAT)/XTOTBV_IN_MESH(JMESH)
00328         ZC_DEPTH_RATIO(JMESH) = ZC_DEPTH_RATIO(JMESH) + XC_DEPTH_RATIO_BV(JCAT)*XBV_IN_MESH(JMESH,JCAT)/XTOTBV_IN_MESH(JMESH)
00329       ENDIF
00330     ENDDO
00331   ENDDO
00332   !
00333   WHERE (ZF_PARAM==0.)
00334     ZF_PARAM=2.5
00335     ZC_DEPTH_RATIO=1.
00336   ENDWHERE
00337   !
00338   !write(*,*) 'f min max isba',MINVAL(ZF_PARAM),MAXVAL(ZF_PARAM)
00339   !write(*,*) 'dc min max isba',MINVAL(ZC_DEPTH_RATIO),MAXVAL(ZC_DEPTH_RATIO)
00340   !
00341   CALL OPEN_FILE('ASCII ',IUNIT,'carte_f_dc.txt','FORMATTED',HACTION='WRITE')
00342   DO JMESH=1,NDIM_FULL
00343     WRITE(IUNIT,*) ZF_PARAM(JMESH),ZC_DEPTH_RATIO(JMESH)
00344   ENDDO
00345   CALL CLOSE_FILE('ASCII ',IUNIT)
00346   !
00347   DEALLOCATE(ZF_PARAM)
00348   DEALLOCATE(ZC_DEPTH_RATIO)
00349   !
00350   WRITE(ILUOUT,*) 'Couplage avec TOPMODEL active'
00351   !
00352 ELSE
00353   !
00354   WRITE(ILUOUT,*) 'Pas de couplage avec TOPMODEL'
00355   !
00356 ENDIF
00357 !
00358 IF (LHOOK) CALL DR_HOOK('PGD_TOPD',1,ZHOOK_HANDLE)
00359 !
00360 END SUBROUTINE PGD_TOPD