treat_field.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 treat_field (UG, U, USS, &
                              HPROGRAM,HSCHEME,HFILETYPE,    &
                              HSUBROUTINE,HFILENAME,HFIELD,   &
                              PPGDARRAY            )  
!     ##############################################################
!
!!**** *TREAT_FIELD* chooses which treatment subroutine to use
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    11/09/95
!!
!!    Modification
!!    25/05/96    (V. Masson) remove useless case for HSUBROUTINE   
!!    29/11/2002  (D. Gazen)  add HSFTYPE argument + call to read_binllvfast routine
!!    03/2004     (V. MAsson) externalization
!!    04/2009     (B. Decharme) Special treatement for gaussian grid
!!    06/2009     (B. Decharme)  call Topographic index statistics calculation
!!    09/2010     (E. Kourzeneva) call reading of the lake database
!!    03/2012     (M. Lafaysse) NETCDF
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
!
USE modd_surf_par, ONLY : xundef
USE modd_pgdwork, ONLY : nsize, nsize_all, xall, nsso_all, xsso_all, &
                         xsumval, xext_all, catype, &
                         xssqo, lssqo, nsso, xmin_work, xmax_work, & 
                         nvalnbr, nvalcount, xvallist, jpvalmax
!
USE modd_surfex_omp, ONLY : nblocktot
USE modd_surfex_mpi, ONLY : nrank, npio, nproc, ncomm, nreq, nindex, idx_r, &
                                nsize_task,nreq, nsize_max=>nsize
!
USE modd_data_lake,      ONLY : ngraddepth_ldb, ngradstatus_ldb 
USE modd_data_cover_par, ONLY : jpcover
!
USE modi_ini_ssowork
USE modi_get_luout
USE modi_read_direct
USE modi_read_direct_gauss
USE modi_read_latlon
USE modi_read_binllv
USE modi_read_binllvfast
USE modi_read_ascllv
USE modi_read_and_send_mpi
USE modi_read_pgd_netcdf
USE modi_average2_mesh
USE modi_make_lcover
USE modi_abor1_sfx
USE modi_average2_cover
USE modi_average2_cti
USE modi_average2_ldb
USE modi_average2_orography
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
 CHARACTER(LEN=6),  INTENT(IN) :: HPROGRAM      ! Type of program
 CHARACTER(LEN=6),  INTENT(IN) :: HSCHEME       ! Scheme treated
 CHARACTER(LEN=6),  INTENT(IN) :: HFILETYPE     ! Type of the data file
 CHARACTER(LEN=6),  INTENT(IN) :: HSUBROUTINE   ! Name of the subroutine to call
 CHARACTER(LEN=28), INTENT(IN) :: HFILENAME     ! Name of the field file.
 CHARACTER(LEN=20), INTENT(IN) :: HFIELD        ! Name of the field.
REAL, DIMENSION(:,:), INTENT(INOUT), OPTIONAL :: PPGDARRAY ! field on MESONH grid
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: I3D_ALL
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZVALLIST, ZVAL
REAL, DIMENSION(:,:), ALLOCATABLE ::  ZEXTVAL
INTEGER, DIMENSION(:,:), ALLOCATABLE :: ISIZE0
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IVALNBR
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: IVALCOUNT
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: ISIZE
!
INTEGER :: IMAX  ! Maximum of times a value has been encountered in the grid mesh
INTEGER :: IVAL  ! Index of this value
#ifdef SFX_MPI
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
INTEGER, DIMENSION(MPI_STATUS_SIZE,NPROC-1) :: ISTATUS2
#endif
INTEGER, DIMENSION(0:NPROC-1) :: ITCOV
INTEGER :: ILUOUT, IS2, INFOMPI, JP, ICPT, JCOV, JI, JL, IREQ, IDX,&
                        IDX_SAVE, JT
LOGICAL :: GMULTITYPE
REAL(KIND=JPRB) :: ZHOOK_HANDLE, ZHOOK_HANDLE_OMP
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TREAT_FIELD_1',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
!*    1.     Selection of type of reading (and point by point treatment)
!            -----------------------------------------------------------
!
gmultitype = .false.
IF (hfiletype=='DIRTYP') gmultitype = .true.
!
SELECT CASE (hfiletype)

  CASE ('DIRECT','DIRTYP')
    IF(ug%G%CGRID=="GAUSS     " .OR. ug%G%CGRID=="IGN       " .OR. ug%G%CGRID=="LONLAT REG")THEN
      CALL read_direct_gauss(ug, u, uss, &
                             hprogram,hscheme,hsubroutine,hfilename,hfield,gmultitype)
    ELSE
      CALL read_direct(ug, u, uss, &
                       hprogram,hscheme,hsubroutine,hfilename,hfield,gmultitype)
    ENDIF

  CASE ('BINLLV')
    IF (hsubroutine=='A_OROG') CALL ini_ssowork
    IF (nrank==npio) CALL read_binllv(ug, u, uss, &
                        hprogram,hsubroutine,hfilename)

  CASE ('BINLLF')
    IF (hsubroutine=='A_OROG') CALL ini_ssowork
    IF (nrank==npio) CALL read_binllvfast(ug, u, uss, &
                            hprogram,hsubroutine,hfilename)

  CASE ('ASCLLV')
    IF (hsubroutine=='A_OROG') CALL ini_ssowork
    IF (nrank==npio) CALL read_ascllv(ug, u, uss, &
                        hprogram,hsubroutine,hfilename)

  CASE ('LATLON')
    IF (hsubroutine=='A_OROG') CALL ini_ssowork
    IF (nrank==npio) CALL read_latlon(ug, u, uss, &
                        hprogram,hscheme,hsubroutine,hfilename)

  CASE ('NETCDF')
    IF (hsubroutine=='A_OROG') CALL ini_ssowork
    IF (nrank==npio) CALL read_pgd_netcdf(ug, u, uss, &
                            hprogram,hscheme,hsubroutine,hfilename,hfield)

  CASE DEFAULT
    CALL abor1_sfx('TREAT_FIELD: FILE TYPE NOT SUPPORTED: '//hfiletype)

END SELECT
!
!-------------------------------------------------------------------------------
!
!nsize contains the number of points found for each of the domain, for each task
ALLOCATE(nsize(u%NSIZE_FULL,SIZE(nsize_all,2)))
!
IF (lhook) CALL dr_hook('TREAT_FIELD_1',1,zhook_handle)
!
IF (nproc>1) THEN
  !
  IF (hfiletype=='DIRECT'.OR.hfiletype=='DIRTYP') THEN
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_21',0,zhook_handle)
    !
    ALLOCATE(isize(nsize_max,SIZE(nsize_all,2),nproc))
    ALLOCATE(isize0(nsize_max,SIZE(nsize_all,2)))
    !
    idx_save = idx_r
    idx = idx_save + nrank
    !
    idx = idx + 1
    isize(:,:,:) = 0
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_21',1,zhook_handle)
    !
!$OMP PARALLEL PRIVATE(ZHOOK_HANDLE_OMP) 
    IF (lhook) CALL dr_hook('TREAT_FIELD_22',0,zhook_handle_omp)
!$OMP DO SCHEDULE(DYNAMIC,1) PRIVATE(JP,ICPT,JI,IREQ,INFOMPI)
    DO jp=0,nproc-1
      IF (jp/=nrank) THEN
        icpt = 0
        DO ji = 1,SIZE(nindex)
          IF (nindex(ji)==jp) THEN
            icpt = icpt + 1
            isize(icpt,:,jp+1) = nsize_all(ji,:)
          ENDIF
        ENDDO
        IF (jp<nrank) THEN
          ireq = jp+1
        ELSE
          ireq = jp
        ENDIF
#ifdef SFX_MPI        
        CALL mpi_isend(isize(:,:,jp+1),SIZE(isize,1)*SIZE(isize,2)*kind(isize)/4,&
                             mpi_integer,jp,idx,ncomm,nreq(ireq),infompi)
#endif
      ENDIF
    ENDDO
!$OMP END DO
    IF (lhook) CALL dr_hook('TREAT_FIELD_22',1,zhook_handle_omp)
!$OMP END PARALLEL
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_23',0,zhook_handle)
    !
    nsize(:,:) = 0
    !for each task
    DO jp=0,nproc-1
      !
      isize0(:,:) = 0
      !
      IF (jp/=nrank) THEN
        !each task receives each ISIZE from each task
#ifdef SFX_MPI        
        CALL mpi_recv(isize0,nsize_max*SIZE(isize0,2)*kind(isize0)/4,mpi_integer,&
                      jp,idx_save+1+jp,ncomm,istatus,infompi)
#endif
      ELSE
        !
        icpt = 0
        DO ji = 1,SIZE(nindex)
          IF (nindex(ji)==jp) THEN
            icpt = icpt + 1
            isize0(icpt,:) = nsize_all(ji,:)
          ENDIF
        ENDDO
        !
      ENDIF
      !
      !nsize is the sum of all parts isize
      nsize(:,:) = nsize(:,:) + isize0(1:nsize_task(nrank),:)
      !
    ENDDO
#ifdef SFX_MPI    
    CALL mpi_waitall(nproc-1,nreq(1:nproc-1),istatus2,infompi)
#endif
    !
    DEALLOCATE(isize,isize0)
    !
    idx_r = idx_r + nproc
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_23',1,zhook_handle)
    !
  ELSE
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_3',0,zhook_handle)
    !
    CALL read_and_send_mpi(nsize_all,nsize)
    !
    IF (lhook) CALL dr_hook('TREAT_FIELD_3',1,zhook_handle)
    !
  ENDIF
  !
ELSE
  IF (lhook) CALL dr_hook('TREAT_FIELD_4',0,zhook_handle)
  nsize(:,:) = nsize_all(:,:)
  IF (lhook) CALL dr_hook('TREAT_FIELD_4',1,zhook_handle)
ENDIF
!
IF (lhook) CALL dr_hook('TREAT_FIELD_5',0,zhook_handle)
!
DEALLOCATE(nsize_all)
!
!
IF (hsubroutine=='A_COVR') THEN
  is2 = SIZE(xall,2)
ELSEIF (hsubroutine=='A_LDBS') THEN
  is2 = ngradstatus_ldb
ELSEIF (hsubroutine=='A_LDBD') THEN
  is2 = ngraddepth_ldb
ELSEIF (hsubroutine=='A_OROG') THEN
  is2 = 2
ELSEIF (hsubroutine=='A_CTI ') THEN
  is2 = 3 
ENDIF
!
!
SELECT CASE (hsubroutine)

  CASE ('A_COVR')
    !
    !to gather the parts of LCOVER sparsed among the tasks
    CALL make_lcover(u%LCOVER)
    !
    !contains the indexes of the covers in XCOVER, associated to their effective
    !num
    ALLOCATE(imask(SIZE(u%LCOVER)))
    imask(:) = 0
    icpt = 0
    DO jcov = 1,SIZE(u%LCOVER)
      IF (u%LCOVER(jcov)) THEN
        icpt = icpt + 1
        imask(jcov) = icpt
      ENDIF
    ENDDO
    !
    !because each task did not necessarily meet the same number of coves, itcov 
    !contains the numbers of covers met for all tasks 
    IF (nproc>1) THEN
#ifdef SFX_MPI
    CALL mpi_allgather(is2,kind(is2)/4,mpi_integer,&
                     itcov,kind(itcov)/4,mpi_integer,ncomm,infompi)
#endif
    ELSE
      itcov(:) = is2
    ENDIF
    !
    !XSUMVAL needs to contain the numbers of times each cover is encountered
    !for the current task 
    ALLOCATE(xsumval(u%NSIZE_FULL,count(u%LCOVER)))
    xsumval(:,:) = 0.
    !
    IF (nproc>1) THEN
      ALLOCATE(zval(u%NSIZE_FULL,maxval(itcov),2))
      DO jp = 0,nproc-1
        !the part of XALL concerning the current task is sent 
        !from each other task
        CALL read_and_send_mpi(xall,zval(:,1:itcov(jp),:),kpio=jp)
        DO jl = 1,itcov(jp)
          DO ji=1,u%NSIZE_FULL
            jcov = nint(zval(ji,jl,1))
            IF (jcov/=0) THEN
              !xsumval2 is the sum of contributions zsumval3 coming from all tasks
              xsumval(ji,imask(jcov)) = xsumval(ji,imask(jcov)) + zval(ji,jl,2)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      DEALLOCATE(zval)
    ELSE
      DO jl = 1,SIZE(xall,2)
        DO ji=1,SIZE(xall,1)
          jcov = nint(xall(ji,jl,1))
          IF (jcov/=0) THEN
            xsumval(ji,imask(jcov)) = xsumval(ji,imask(jcov)) + xall(ji,jl,2)
          ENDIF
        ENDDO
      ENDDO
    ENDIF
    DEALLOCATE(xall,imask)
    !
    !
  CASE ('A_LDBD','A_LDBS','A_OROG','A_CTI ')
    !
    !XSUMVAL needs to contain the numbers of times each quantity is encountered
    !for the current task    
    ALLOCATE(xsumval(u%NSIZE_FULL,is2))
    xsumval(:,:) = 0.
    IF (nproc>1) THEN
      ALLOCATE(zval(u%NSIZE_FULL,is2,1))
      DO jp = 0,nproc-1
        CALL read_and_send_mpi(xall(:,:,1),zval(:,:,1),kpio=jp)
        !sum of contributions coming from all tasks
        xsumval(:,:) = xsumval(:,:) + zval(:,:,1)
      ENDDO
      DEALLOCATE(zval)
    ELSE
      xsumval(:,:) = xall(:,:,1)
    ENDIF
    DEALLOCATE(xall)
    !
    !
    IF (hsubroutine=='A_OROG' .OR. hsubroutine=='A_CTI ') THEN
      !
      !special fields
      !
      IF (hsubroutine=='A_CTI ') THEN
        !max and min
        IF (nproc>1) THEN
          ALLOCATE(zextval(u%NSIZE_FULL,1))
          DO jp = 0,nproc-1
            CALL read_and_send_mpi(xext_all(:,1),zextval(:,1),kpio=jp)
            xmax_work(:) = max(xmax_work,zextval(:,1))
          ENDDO
          DO jp = 0,nproc-1
            CALL read_and_send_mpi(xext_all(:,2),zextval(:,1),kpio=jp)
            xmin_work(:) = min(xmin_work,zextval(:,1))
          ENDDO   
          DEALLOCATE(zextval)
        ELSE
          xmax_work(:) = xext_all(:,1)
          xmin_work(:) = xext_all(:,2)
        ENDIF
        !
      ELSEIF (hsubroutine=='A_OROG') THEN
        !max and min
        IF (nproc>1) THEN
          ALLOCATE(zextval(u%NSIZE_FULL,1))
          DO jp = 0,nproc-1
            CALL read_and_send_mpi(xext_all(:,1),zextval(:,1),kpio=jp)
            uss%XMAX_ZS(:) = max(uss%XMAX_ZS,zextval(:,1))
          ENDDO 
          DO jp = 0,nproc-1
            CALL read_and_send_mpi(xext_all(:,2),zextval(:,1),kpio=jp)
            uss%XMIN_ZS(:) = min(uss%XMIN_ZS,zextval(:,1))
          ENDDO   
          DEALLOCATE(zextval)
        ELSE
          uss%XMAX_ZS(:) = xext_all(:,1)
          uss%XMIN_ZS(:) = xext_all(:,2)
        ENDIF
        !
        !sso fields   
        ALLOCATE(xssqo(u%NSIZE_FULL,nsso,nsso))
        xssqo(:,:,:) = -xundef
        IF (nproc>1) THEN
          ALLOCATE(zval(u%NSIZE_FULL,nsso,nsso))
          DO jp = 0,nproc-1
            CALL read_and_send_mpi(xsso_all,zval,kpio=jp)
            xssqo(:,:,:) = max(xssqo(:,:,:),zval)
          ENDDO
          DEALLOCATE(zval)
        ELSE
          xssqo(:,:,:) = xsso_all(:,:,:)
        ENDIF
        DEALLOCATE(xsso_all)
        !
        ALLOCATE(lssqo(u%NSIZE_FULL,nsso,nsso))
        lssqo(:,:,:) = .false.
        IF (nproc>1) THEN
          ALLOCATE(i3d_all(u%NSIZE_FULL,nsso,nsso)) 
          DO jp = 0,nproc-1       
            CALL read_and_send_mpi(nsso_all,i3d_all,kpio=jp)
            WHERE (i3d_all(:,:,:)==1) lssqo(:,:,:) = .true.
          ENDDO
          DEALLOCATE(i3d_all)
        ELSE
          WHERE(nsso_all(:,:,:)==1) lssqo(:,:,:) = .true.
        ENDIF
        DEALLOCATE(nsso_all)
      ENDIF
      DEALLOCATE(xext_all)
      !
    ENDIF
    !
    !
  CASE ('A_MESH')
   IF (catype/='MAJ') THEN

    ALLOCATE(xsumval(u%NSIZE_FULL,SIZE(xall,2)))
    !most simple case
    IF (nproc>1) THEN
      xsumval(:,:) = 0.
      ALLOCATE(zval(u%NSIZE_FULL,SIZE(xall,2),1))
      DO jp = 0,nproc-1
        CALL read_and_send_mpi(xall(:,:,1),zval(:,:,1),kpio=jp)
        xsumval(:,:) = xsumval(:,:) + zval(:,:,1)
      ENDDO
      DEALLOCATE(zval)
    ELSE
      xsumval(:,:) = xall(:,:,1)
    ENDIF
    DEALLOCATE(xall)

   ELSE

     ALLOCATE(xsumval(u%NSIZE_FULL,1))
     IF (hfiletype=='DIRECT' .AND. nproc>1) THEN
       CALL abor1_sfx("TREAT_FIELD: MAJ is not possible with DIRECT filetype and NPROC>1")
     ELSE
       ALLOCATE(ivalnbr(u%NSIZE_FULL,SIZE(nvalnbr,2)),ivalcount(u%NSIZE_FULL,jpvalmax,SIZE(nvalnbr,2)),&
                zvallist(u%NSIZE_FULL,jpvalmax,SIZE(nvalnbr,2)))
       IF (nproc>1) THEN
         CALL read_and_send_mpi(nvalnbr,ivalnbr)
         CALL read_and_send_mpi(nvalcount,ivalcount)
         CALL read_and_send_mpi(xvallist,zvallist)
       ELSE
         ivalnbr = nvalnbr
         ivalcount = nvalcount
         zvallist = xvallist
       ENDIF
       DEALLOCATE(nvalnbr,nvalcount,xvallist)
       DO jt=1,SIZE(nsize,2)
         DO ji=1,SIZE(nsize,1)
           IF(nsize(ji,jt)==0) cycle
             !* determines the index of the value which has been the most encountered
             !  in the grid mesh
             imax=0
             DO jl=1,ivalnbr(ji,jt)
               IF (ivalcount(ji,jl,jt)>imax) THEN
                 imax = ivalcount(ji,jl,jt)
                 ival = jl
               END IF
             END DO
             !* sets this value to the PGD field
             xsumval(ji,jt)=zvallist(ji,ival,jt)
         END DO
       ENDDO
       DEALLOCATE(ivalnbr,ivalcount,zvallist)
     ENDIF
   ENDIF
    !
END SELECT
!
IF (lhook) CALL dr_hook('TREAT_FIELD_5',1,zhook_handle)
!-------------------------------------------------------------------------------
!
!*    2.     Call to the adequate subroutine (global treatment)
!            --------------------------------------------------
!
IF (lhook) CALL dr_hook('TREAT_FIELD_6',0,zhook_handle)
!
SELECT CASE (hsubroutine)

  CASE ('A_COVR')
    CALL average2_cover(u, hprogram)

  CASE ('A_OROG')
    CALL average2_orography(uss)

  CASE ('A_CTI ')
    CALL average2_cti

  CASE ('A_LDBD')
    CALL average2_ldb(ppgdarray(:,1),'D',1)

  CASE ('A_LDBS')
    CALL average2_ldb(ppgdarray(:,1),'S',1)
    
  CASE ('A_MESH')
    IF (.NOT. PRESENT(ppgdarray)) THEN
      WRITE(iluout,*) 'You asked to average a PGD field with A_MESH option,'
      WRITE(iluout,*) 'but you did not give the array to store this field'
      CALL abor1_sfx('TREAT_FIELD: ARRAY IS MISSING')
    END IF
    CALL average2_mesh(ppgdarray)

END SELECT
!
IF (lhook) CALL dr_hook('TREAT_FIELD_6',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE treat_field