From e7ed4812e85105c38b86d9f62a385a07d0d7f093 Mon Sep 17 00:00:00 2001 From: Juan Escobar <juan.escobar@aero.obs-mip.fr> Date: Tue, 24 Nov 2015 09:36:09 +0000 Subject: [PATCH] Juan 24/11/2015: add for PREPLL --- src/LIB/SURCOUCHE/src/update_nhalo1d.f90 | 384 ++++++++++++++++++ src/MNH/extend_grid_parameter_mnh.f90 | 176 ++++++++ src/MNH/goto_model_surfex_mnh.f90 | 168 ++++++++ src/MNH/mnh_surf_grid_io_init.f90 | 158 +++++++ src/SURFEX/extend_grid_on_halo.F90 | 87 ++++ src/SURFEX/extend_grid_on_halo_cartesian.F90 | 117 ++++++ src/SURFEX/extend_grid_on_halo_conf_proj.F90 | 117 ++++++ src/SURFEX/goto_model_mnh.F90 | 141 +++++++ src/SURFEX/hor_interpol_1cov.F90 | 133 ++++++ src/SURFEX/hor_interpol_arome_1cov.F90 | 105 +++++ src/SURFEX/hor_interpol_buffer_1cov.F90 | 78 ++++ src/SURFEX/hor_interpol_cartesian_1cov.F90 | 107 +++++ src/SURFEX/hor_interpol_conf_proj_1cov.F90 | 142 +++++++ src/SURFEX/hor_interpol_gauss_1cov.F90 | 254 ++++++++++++ src/SURFEX/hor_interpol_latlon_1cov.F90 | 91 +++++ src/SURFEX/hor_interpol_none_1cov.F90 | 78 ++++ src/SURFEX/hor_interpol_rotlatlon_1cov.F90 | 305 ++++++++++++++ src/SURFEX/mode_extend_grid_parameter.F90 | 117 ++++++ .../read_covers_and_av_pgd_on_layers.F90 | 210 ++++++++++ src/SURFEX/write_lcover.F90 | 86 ++++ 20 files changed, 3054 insertions(+) create mode 100644 src/LIB/SURCOUCHE/src/update_nhalo1d.f90 create mode 100644 src/MNH/extend_grid_parameter_mnh.f90 create mode 100644 src/MNH/goto_model_surfex_mnh.f90 create mode 100644 src/MNH/mnh_surf_grid_io_init.f90 create mode 100644 src/SURFEX/extend_grid_on_halo.F90 create mode 100644 src/SURFEX/extend_grid_on_halo_cartesian.F90 create mode 100644 src/SURFEX/extend_grid_on_halo_conf_proj.F90 create mode 100644 src/SURFEX/goto_model_mnh.F90 create mode 100644 src/SURFEX/hor_interpol_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_arome_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_buffer_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_cartesian_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_conf_proj_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_gauss_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_latlon_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_none_1cov.F90 create mode 100644 src/SURFEX/hor_interpol_rotlatlon_1cov.F90 create mode 100644 src/SURFEX/mode_extend_grid_parameter.F90 create mode 100644 src/SURFEX/read_covers_and_av_pgd_on_layers.F90 create mode 100644 src/SURFEX/write_lcover.F90 diff --git a/src/LIB/SURCOUCHE/src/update_nhalo1d.f90 b/src/LIB/SURCOUCHE/src/update_nhalo1d.f90 new file mode 100644 index 000000000..48c3566b6 --- /dev/null +++ b/src/LIB/SURCOUCHE/src/update_nhalo1d.f90 @@ -0,0 +1,384 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ################################################################ + SUBROUTINE UPDATE_NHALO1D( NHALO, PFIELD1D, KISIZE_ll, KJSIZE_ll, KXOR, KXEND, KYOR, KYEND, HREC ) +! ################################################################ +! +!!**** *UPDATE_NHALO1D* - routine to update a halo of size NHALO for 1D field PFIELD1D +!! of physical size KSIZE +!! WARNING : if NHALO > KISIZE_ll or NHALO > KJSIZE_ll, this routine will crash +!! +!! PURPOSE +!! ------- +!! +!!** METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge *LA - CNRS* +!! +!! MODIFICATIONS +!! ------------- +!! Original 05/2015 +!! M.Moge 07/2015 bug fix : TAG computation and ICARDDIF computation +!! + use a temp field PFIELDTMP for SEND_RECV_FIELD +!! M.Moge 08/2015 calling ABORT if local subdomain is of size < NHALO +!! (this causes problems on the boundary of the domain) +!! M.Moge 08/2015 bug fix : changing the computation of IISIZE +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_SURF_PAR, ONLY : NUNDEF +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +USE MODE_ll +USE MODE_EXCHANGE_ll, ONLY : SEND_RECV_FIELD +USE MODE_SPLITTING_ll, ONLY : SPLIT2 +USE MODD_VAR_ll, ONLY : NPROC, IP, YSPLITTING, NMNH_COMM_WORLD +USE MODD_MPIF +USE MODD_STRUCTURE_ll, ONLY : ZONE_ll, CRSPD_ll +USE MODE_TOOLS_ll, ONLY : INTERSECTION +USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments +! ------------------------- +! +INTEGER, INTENT(IN) :: NHALO ! size of HALO to update +REAL, DIMENSION(:), INTENT(INOUT) :: PFIELD1D ! field to be updated +INTEGER, INTENT(IN) :: KISIZE_ll ! global physical size of the domain in X direction +INTEGER, INTENT(IN) :: KJSIZE_ll ! global physical size of the domain in Y direction +INTEGER, INTENT(IN) :: KXOR ! origin and +INTEGER, INTENT(IN) :: KXEND ! end of local subdomain +INTEGER, INTENT(IN) :: KYOR ! in global domain +INTEGER, INTENT(IN) :: KYEND ! defined by KISIZE_ll and KJSIZE_ll +CHARACTER(LEN=6), INTENT(IN) :: HREC ! name of the parameter, 'XX' or 'YY' +! +!* 0.2 Declarations of local variables +! ------------------------------- +! +REAL, DIMENSION(:,:,:), ALLOCATABLE :: PFIELDTMP ! PFIELDTMP will get the communicated halos +! +!* other variables +! +INTEGER :: JI,JJ ! loop controls +REAL(KIND=JPRB) :: ZHOOK_HANDLE +INTEGER :: IINFO_ll +! +! structures for the partitionning +! +TYPE(ZONE_ll), DIMENSION(NPROC) :: TZSPLITTING_PHYS !physical splitting of the field +TYPE(ZONE_ll), DIMENSION(NPROC) :: TZSPLITTING_EXT !extended splitting of the field +TYPE(ZONE_ll) :: TZFIELD_ll ! global field +! +! structures for the communications +! +TYPE(ZONE_ll), ALLOCATABLE, DIMENSION(:) :: TZSEND, TZRECV +TYPE(CRSPD_ll), POINTER :: TZCRSPDSEND, TZCRSPDRECV +TYPE(CRSPD_ll), ALLOCATABLE, DIMENSION(:), TARGET :: TZCRSPDSENDTAB, TZCRSPDRECVTAB +! +INTEGER :: J +INTEGER :: INBMSG +INTEGER :: ICARD +INTEGER :: ICARDDIF +INTEGER :: IISIZE, IJSIZE, IKSIZE ! local sizes of the field +INTEGER :: IXOR_ll, IYOR_ll, IZOR_ll, IXEND_ll, IYEND_ll, IZEND_ll ! origin and end of local physical subdomain +REAL, DIMENSION(:,:,:), ALLOCATABLE :: PFIELD3D ! field to be updated +INTEGER , DIMENSION(NPROC) :: IXORARRAY_ALL, IYORARRAY_ALL, IXENDARRAY_ALL, IYENDARRAY_ALL ! array containing origin and end of each subdomain in local coordinates +INTEGER , DIMENSION(NPROC) :: IDISPLS !displacement array for MPI_ALLGATHERV +INTEGER , DIMENSION(NPROC) :: IRECVCOUNTS !nteger array containing the number of elements that are to be received from each process +!------------------------------------------------------------------------------ +! +!* 1. Coherence check +! -------------------------------------------------------------- +! +! verify that the local sizes of PFIELD3D, the halo size (NHALO) +! and the global physical sizes (KISIZE_ll,KJSIZE_ll) are OK +!IF ( ) THEN +! +! RETURN +!ENDIF +!! +!! we assume that the sizes are correct +!! +!------------------------------------------------------------------------------ +! +!* 1. Partitionning of the field +! -------------------------------------------------------------- +! +! Si le sous-domaine local est plus petit que NHALO, alors on aura des problemes sur les sous-domaines +! voisins des sous-domaines au bord : +! si le sous-domaine local est a moins de NHALO points du bord mais n'est pas au bord, +! alors il y a un probleme car on ne recevra pas de valeurs pour les points du HALO situes sur le halo externe +! Donc on fait un WARNING et un ABORT +! +IF ( NHALO > KXEND - KXOR + 1 .OR. NHALO > KYEND - KYOR + 1 ) THEN + WRITE(*,*) "ERROR in UPDATE_NHALO1D : size of local subdomain is (", KXEND - KXOR + 1,",",KYEND - KYOR + 1, \ + ") which is less than NHALO=",NHALO + WRITE(*,*) "Try with less MPI processes or a larger domain" + CALL ABORT +ENDIF +! +! physical splitting of the field +! +IXOR_ll = KXOR + NHALO - JPHEXT +IXEND_ll = KXEND + NHALO - JPHEXT +IYOR_ll = KYOR + NHALO - JPHEXT +IYEND_ll = KYEND + NHALO - JPHEXT +IZOR_ll = 1 +IZEND_ll = 1 +! +IISIZE = IXEND_ll - IXOR_ll + 1 + 2 * NHALO +IJSIZE = IYEND_ll - IYOR_ll + 1 + 2 * NHALO +IKSIZE = 1 +ALLOCATE( PFIELDTMP(IISIZE,IJSIZE,IKSIZE) ) +ALLOCATE( PFIELD3D(IISIZE,IJSIZE,IKSIZE) ) +IF ( HREC == 'XX' ) THEN + DO JI = 1, IJSIZE + PFIELD3D(:,JI,1) = PFIELD1D(:) + ENDDO +ELSE IF ( HREC == 'YY' ) THEN + DO JI = 1, IISIZE + PFIELD3D(JI,:,1) = PFIELD1D(:) + ENDDO +ENDIF +! +! we get the splitting of the physical domain with a MPI_AllGatherv +! +DO JI = 1, NPROC + IDISPLS( JI ) = JI-1 +ENDDO +IRECVCOUNTS(:) = 1 +CALL MPI_ALLGATHERV( KXOR, 1, MPI_INTEGER, IXORARRAY_ALL, IRECVCOUNTS, IDISPLS, MPI_INTEGER, NMNH_COMM_WORLD, IINFO_ll) +CALL MPI_ALLGATHERV( KXEND, 1, MPI_INTEGER, IXENDARRAY_ALL, IRECVCOUNTS, IDISPLS, MPI_INTEGER, NMNH_COMM_WORLD, IINFO_ll) +CALL MPI_ALLGATHERV( KYOR, 1, MPI_INTEGER, IYORARRAY_ALL, IRECVCOUNTS, IDISPLS, MPI_INTEGER, NMNH_COMM_WORLD, IINFO_ll) +CALL MPI_ALLGATHERV( KYEND, 1, MPI_INTEGER, IYENDARRAY_ALL, IRECVCOUNTS, IDISPLS, MPI_INTEGER, NMNH_COMM_WORLD, IINFO_ll) +! +DO JI = 1, NPROC + TZSPLITTING_PHYS(JI)%NUMBER = JI + TZSPLITTING_PHYS(JI)%NXOR = IXORARRAY_ALL(JI) + NHALO - JPHEXT + TZSPLITTING_PHYS(JI)%NXEND = IXENDARRAY_ALL(JI) + NHALO - JPHEXT + TZSPLITTING_PHYS(JI)%NYOR = IYORARRAY_ALL(JI) + NHALO - JPHEXT + TZSPLITTING_PHYS(JI)%NYEND = IYENDARRAY_ALL(JI) + NHALO - JPHEXT + TZSPLITTING_PHYS(JI)%NZOR = 1 + TZSPLITTING_PHYS(JI)%NZEND = 1 +ENDDO +! +! extended splitting of the field +! +DO JI = 1, NPROC + TZSPLITTING_EXT(JI)%NUMBER = JI + TZSPLITTING_EXT(JI)%NXOR = TZSPLITTING_PHYS(JI)%NXOR - NHALO + TZSPLITTING_EXT(JI)%NXEND = TZSPLITTING_PHYS(JI)%NXEND + NHALO + TZSPLITTING_EXT(JI)%NYOR = TZSPLITTING_PHYS(JI)%NYOR - NHALO + TZSPLITTING_EXT(JI)%NYEND = TZSPLITTING_PHYS(JI)%NYEND + NHALO + TZSPLITTING_EXT(JI)%NZOR = 1 + TZSPLITTING_EXT(JI)%NZEND = 1 +ENDDO +! +!------------------------------------------------------------------------------ +! +!* 2. Creation of global ZONE_ll for the whole field +! -------------------------------------------------------------- +! + TZFIELD_ll%NXOR = 1 + TZFIELD_ll%NYOR = 1 + TZFIELD_ll%NZOR = 1 + TZFIELD_ll%NXEND = KISIZE_ll + 2 * NHALO + IF( KJSIZE_ll > 1 ) THEN + TZFIELD_ll%NYEND = KJSIZE_ll + 2 * NHALO + ELSE + TZFIELD_ll%NYEND = KJSIZE_ll + ENDIF + TZFIELD_ll%NZEND = 1 +!------------------------------------------------------------------------------ +! +!* 3. Preparing the structures for the communications +! -------------------------------------------------------------- +! + ! + ! ######## initializing the structures for the SEND ######## + ! we take the intersection of the local physical subdomain with all extended subdomains + ! + ALLOCATE(TZSEND(NPROC)) + CALL INTERSECTION( TZSPLITTING_EXT, NPROC, TZSPLITTING_PHYS(IP), TZSEND) + ! il faut initialiser le TAG de manière a avoir un meme tag unique pour le send et le recv : + ! on concatene le num du proc qui envoie et le num du proc qui recoit + DO J = 1, NPROC + IF ( TZSEND(J)%NUMBER > 0 ) THEN + TZSEND(J)%MSSGTAG = IP * 10**(FLOOR(LOG10(real(TZSEND(J)%NUMBER)))+1) + TZSEND(J)%NUMBER + ENDIF + ENDDO + ! switching to local coordinates + DO J = 1, NPROC + IF ( TZSEND(J)%NUMBER > 0 ) THEN + TZSEND(J)%NXOR = TZSEND(J)%NXOR - IXOR_ll + NHALO + 1 + TZSEND(J)%NXEND = TZSEND(J)%NXEND - IXOR_ll + NHALO + 1 + TZSEND(J)%NYOR = TZSEND(J)%NYOR - IYOR_ll + 1 + TZSEND(J)%NYEND = TZSEND(J)%NYEND - IYOR_ll + 1 + IF( KJSIZE_ll > 1 ) THEN + TZSEND(J)%NYOR = TZSEND(J)%NYOR + NHALO + TZSEND(J)%NYEND = TZSEND(J)%NYEND + NHALO + ENDIF + TZSEND(J)%NZOR = 1 + TZSEND(J)%NZEND = 1 + ENDIF + ENDDO + ! we do not need the Z dimension + DO J = 1, NPROC + IF ( TZSEND(J)%NUMBER > 0 ) THEN + TZSEND(J)%NZOR = 1 + TZSEND(J)%NZEND = 1 + ENDIF + ENDDO + ! switching from an array of CRSPD_ll to a CRSPD_ll pointer + INBMSG = 0 + DO J = 1, NPROC + IF ( TZSEND(J)%NUMBER > 0 ) THEN + INBMSG = INBMSG+1 + ENDIF + ENDDO + ALLOCATE( TZCRSPDSENDTAB(INBMSG) ) + ICARD = 0 + ICARDDIF = 0 + DO J = 1, NPROC + IF ( TZSEND(J)%NUMBER > 0 ) THEN + ICARD = ICARD+1 + IF ( TZSEND(J)%NUMBER /= IP ) THEN + ICARDDIF = ICARDDIF+1 + ENDIF + TZCRSPDSENDTAB(ICARD)%TELT = TZSEND(J) + IF ( ICARD == INBMSG ) THEN + TZCRSPDSENDTAB(ICARD)%TNEXT => NULL() + ELSE + TZCRSPDSENDTAB(ICARD)%TNEXT => TZCRSPDSENDTAB(ICARD+1) + ENDIF + ENDIF + ENDDO + DO J = 1, ICARD + TZCRSPDSENDTAB(J)%NCARD = ICARD-J+1 + IF ( TZCRSPDSENDTAB(J)%TELT%NUMBER > IP ) THEN + TZCRSPDSENDTAB(J)%NCARDDIF = ICARDDIF-J+2 + ELSE + TZCRSPDSENDTAB(J)%NCARDDIF = ICARDDIF-J+1 + ENDIF + ENDDO + IF (ICARD > 0) THEN + TZCRSPDSEND => TZCRSPDSENDTAB(1) + ELSE + TZCRSPDSEND => NULL() + ENDIF + ! + ! ######## initializing the structures for the RECV ######## + ! we take the intersection of the local extend subdomain with all physical subdomains + ! + ALLOCATE(TZRECV(NPROC)) + CALL INTERSECTION( TZSPLITTING_PHYS, NPROC, TZSPLITTING_EXT(IP), TZRECV) + ! il faut initialiser le TAG de manière a avoir un meme tag unique pour le send et le recv : + ! on concatene le num du proc qui envoie et le num du proc qui recoit + DO J = 1, NPROC + IF ( TZRECV(J)%NUMBER > 0 ) THEN + TZRECV(J)%MSSGTAG = TZRECV(J)%NUMBER * 10**(FLOOR(LOG10(real(IP)))+1) + IP + ENDIF + ENDDO + ! switching to local coordinates + DO J = 1, NPROC + IF ( TZRECV(J)%NUMBER > 0 ) THEN + TZRECV(J)%NXOR = TZRECV(J)%NXOR - IXOR_ll + NHALO + 1 + TZRECV(J)%NXEND = TZRECV(J)%NXEND - IXOR_ll + NHALO + 1 + TZRECV(J)%NYOR = TZRECV(J)%NYOR - IYOR_ll + 1 + TZRECV(J)%NYEND = TZRECV(J)%NYEND - IYOR_ll + 1 + IF( KJSIZE_ll > 1 ) THEN + TZRECV(J)%NYOR = TZRECV(J)%NYOR + NHALO + TZRECV(J)%NYEND = TZRECV(J)%NYEND + NHALO + ENDIF + TZRECV(J)%NZOR = 1 + TZRECV(J)%NZEND = 1 + ENDIF + ENDDO + ! we do not need the Z dimension + DO J = 1, NPROC + IF ( TZRECV(J)%NUMBER > 0 ) THEN + TZRECV(J)%NZOR = 1 + TZRECV(J)%NZEND = 1 + ENDIF + ENDDO + ! switching from an array of CRSPD_ll to a CRSPD_ll pointer + INBMSG = 0 + DO J = 1, NPROC + IF ( TZRECV(J)%NUMBER > 0 ) THEN + INBMSG = INBMSG+1 + ENDIF + ENDDO + ALLOCATE( TZCRSPDRECVTAB(INBMSG) ) + ICARD = 0 + ICARDDIF = 0 + DO J = 1, NPROC + IF ( TZRECV(J)%NUMBER > 0 ) THEN + ICARD = ICARD+1 + IF ( TZRECV(J)%NUMBER /= IP ) THEN + ICARDDIF = ICARDDIF+1 + ENDIF + TZCRSPDRECVTAB(ICARD)%TELT = TZRECV(J) + IF ( ICARD == INBMSG ) THEN + TZCRSPDRECVTAB(ICARD)%TNEXT => NULL() + ELSE + TZCRSPDRECVTAB(ICARD)%TNEXT => TZCRSPDRECVTAB(ICARD+1) + ENDIF + ENDIF + ENDDO + DO J = 1, ICARD + TZCRSPDRECVTAB(J)%NCARD = ICARD-J+1 + IF ( TZCRSPDRECVTAB(J)%TELT%NUMBER > IP ) THEN + TZCRSPDRECVTAB(J)%NCARDDIF = ICARDDIF-J+2 + ELSE + TZCRSPDRECVTAB(J)%NCARDDIF = ICARDDIF-J+1 + ENDIF + ENDDO + IF (ICARD > 0) THEN + TZCRSPDRECV => TZCRSPDRECVTAB(1) + ELSE + TZCRSPDRECV => NULL() + ENDIF +! +!------------------------------------------------------------------------------ +! +!* 5. Communications +! ------------------------------------------------------------ +! +PFIELDTMP(:,:,:) = PFIELD3D(:,:,:) +CALL SEND_RECV_FIELD( TZCRSPDSEND, TZCRSPDRECV, PFIELD3D, PFIELDTMP, IINFO_ll) +! +IF ( HREC == 'XX' ) THEN + PFIELD1D(:) = PFIELDTMP(:,NHALO+1,1) +ELSE IF ( HREC == 'YY' ) THEN + PFIELD1D(:) = PFIELDTMP(NHALO+1,:,1) +ENDIF +DEALLOCATE(PFIELD3D) +DEALLOCATE(PFIELDTMP) +DEALLOCATE(TZSEND) +DEALLOCATE(TZRECV) +IF (LHOOK) CALL DR_HOOK('UPDATE_NHALO1D',1,ZHOOK_HANDLE) +!--------------------------------------------------------------------------- +! +END SUBROUTINE UPDATE_NHALO1D diff --git a/src/MNH/extend_grid_parameter_mnh.f90 b/src/MNH/extend_grid_parameter_mnh.f90 new file mode 100644 index 000000000..9a442818e --- /dev/null +++ b/src/MNH/extend_grid_parameter_mnh.f90 @@ -0,0 +1,176 @@ +! ############################################################# + SUBROUTINE EXTEND_GRID_PARAMETERX1_MNH(HGRID,HREC,KDIM,KSIZE,KIMAX,KJMAX,PFIELD,PFIELD_EXTEND) +! ############################################################# +! +!!**** * - routine to extend a real splitted array on SURFEX halo +! +! Author +! M.Moge 01/03/2015 +! +! Modifications +! 06/2015 (M.Moge) using MPI_BCAST to have the same space step on all processes + calling UPDATE_NHALO1D +! 07/2015 (M.Moge) initializing ZY and ZY to zero +! 08/2015 (M.Moge) bug fix in the call to UPDATE_NHALO1D : IIMAX_ll instead of IJMAX_ll +! +USE MODD_IO_SURF_MNH, ONLY : NHALO +USE MODD_VAR_ll, ONLY : NPROC, IP, MPI_PRECISION, NMNH_COMM_WORLD +USE MODD_MPIF +USE MODE_TOOLS_ll, ONLY : INTERSECTION, LWEST_ll, LEAST_ll, LNORTH_ll, LSOUTH_ll +USE MODI_UPDATE_NHALO1D +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments +! +CHARACTER(LEN=10), INTENT(IN) :: HGRID ! grid type +CHARACTER(LEN=6), INTENT(IN) :: HREC ! name of the parameter +INTEGER, INTENT(IN) :: KDIM ! size of PFIELD +INTEGER, INTENT(IN) :: KSIZE ! size of PFIELD_EXTEND +INTEGER, INTENT(IN) :: KIMAX !(local) dimension of the domain - X direction +INTEGER, INTENT(IN) :: KJMAX !(local) dimension of the domain - Y direction +REAL, DIMENSION(KDIM ), INTENT(IN) :: PFIELD ! real field for complete grid +REAL, DIMENSION(KSIZE), INTENT(OUT):: PFIELD_EXTEND! real field for splitted grid +! +!* 0.2 Declarations of local variables +! +INTEGER :: JI, JJ +INTEGER :: IIB, IIE, IJB, IJE +INTEGER :: NIMAX, NJMAX !(local) dimensions of the subdomain +INTEGER :: IXOR, IYOR !origin of local subdomain +! +REAL, DIMENSION(:), ALLOCATABLE :: ZX, ZY +REAL :: ZDX, ZDY +! +INTEGER :: IINDEX +INTEGER :: IINFO_ll +INTEGER :: IIMAX_ll, IJMAX_ll +!------------------------------------------------------------------------------- +! +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +CALL GET_OR_ll ('B',IXOR,IYOR) +NIMAX=IIE-IIB+1 +NJMAX=IJE-IJB+1 +! +CALL GET_GLOBALDIMS_ll(IIMAX_ll, IJMAX_ll) +! +IF (HREC=='XX' .OR. HREC=='DX') THEN + ALLOCATE(ZX(NIMAX+2*NHALO)) + ZX(:) = 0. + ZX(1+NHALO:NIMAX+NHALO) = PFIELD(1:NIMAX) + IF (HREC=='DX') THEN + ZDX = PFIELD(1) + DO JI=NHALO,1,-1 + ZX(JI) = ZDX + ZX(NIMAX+2*NHALO-JI+1) = ZDX + END DO + ELSE IF (HREC=='XX') THEN + ! UPDATE_NHALO1D should be called after the extrapolation on the outer boundaries + ! to avoid any problem when the halo is larger than the neighbouring subdomain + ! but this case is not handled and will lead to an ABORT in UPDATE_NHALO1D + CALL UPDATE_NHALO1D( NHALO, ZX, IIMAX_ll, IJMAX_ll, IIB+IXOR-1,IIE+IXOR-1,IJB+IYOR-1,IJE+IYOR-1, HREC) + ZDX = 0. + IF (NIMAX>1) ZDX = PFIELD(2) - PFIELD(1) + IF (NIMAX==1) ZDX = PFIELD(1) ! in 1D conf, one assumes that grid + ! is located between X=DX/2 and X=3DX/2 + CALL MPI_BCAST(ZDX, 1, MPI_PRECISION, 0, NMNH_COMM_WORLD, IINFO_ll) + IF( LWEST_ll() ) THEN + DO JI=NHALO,1,-1 + ZX(JI) = ZX(JI+1) - ZDX + END DO + ENDIF + IF( LEAST_ll() ) THEN + DO JI=NHALO,1,-1 + ZX(NIMAX+2*NHALO-JI+1) = ZX(NIMAX+2*NHALO-JI) + ZDX + END DO + ENDIF + END IF + +! + DO JJ=1,NJMAX+2*NHALO + DO JI=1,NIMAX+2*NHALO + IINDEX = JI+(JJ-1)*(NIMAX+2*NHALO) + PFIELD_EXTEND(IINDEX) = ZX(JI) + END DO + END DO + DEALLOCATE(ZX) + +ELSEIF (HREC=='YY' .OR. HREC=='DY') THEN + ALLOCATE(ZY(NJMAX+2*NHALO)) + ZY(:) = 0. + DO JJ=1+NHALO,NJMAX+NHALO + ZY(JJ) = PFIELD(1+(JJ-NHALO-1)*KIMAX) + END DO + IF (HREC=='DY') THEN + ZDY = PFIELD(1) + DO JJ=NHALO,1,-1 + ZY(JJ) = ZDY + ZY(NJMAX+2*NHALO-JJ+1) = ZDY + END DO + ELSE IF (HREC=='YY') THEN + ! UPDATE_NHALO1D should be called after the extrapolation on the outer boundaries + ! to avoid any problem when the halo is larger than the neighbouring subdomain + ! but this case is not handled and will lead to an ABORT in UPDATE_NHALO1D + CALL UPDATE_NHALO1D( NHALO, ZY, IIMAX_ll, IJMAX_ll, IIB+IXOR-1,IIE+IXOR-1,IJB+IYOR-1,IJE+IYOR-1, HREC) + ZDY = 0. + IF (NJMAX>1) ZDY = PFIELD(1+KIMAX) - PFIELD(1) + IF (NJMAX==1) ZDY = PFIELD(1) ! in 1D or 2D conf, one assumes that grid + ! is located between Y=DY/2 and Y=3DY/2 + CALL MPI_BCAST(ZDY, 1, MPI_PRECISION, 0, NMNH_COMM_WORLD, IINFO_ll) + IF( LSOUTH_ll() ) THEN + DO JJ=NHALO,1,-1 + ZY(JJ) = ZY(JJ+1) - ZDY + END DO + ENDIF + IF( LNORTH_ll() ) THEN + DO JJ=NHALO,1,-1 + ZY(NJMAX+2*NHALO-JJ+1) = ZY(NJMAX+2*NHALO-JJ) + ZDY + END DO + ENDIF + END IF + + DO JJ=1,NJMAX+2*NHALO + DO JI=1,NIMAX+2*NHALO + IINDEX = JI+(JJ-1)*(NIMAX+2*NHALO) + PFIELD_EXTEND(IINDEX) = ZY(JJ) + END DO + END DO + DEALLOCATE(ZY) + +END IF +! +!------------------------------------------------------------------------------- +END SUBROUTINE EXTEND_GRID_PARAMETERX1_MNH +! +! +! ############################################################# + SUBROUTINE EXTEND_GRID_PARAMETERN0_MNH(HGRID,HREC,KFIELD,KFIELD_EXTEND) +! ############################################################# +! +!!**** * - routine to "extend" an integer related to splitted grid on SURFEX halo +! +! +! +USE MODE_ll +! +USE MODD_IO_SURF_MNH, ONLY : NHALO +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments +! +CHARACTER(LEN=10), INTENT(IN) :: HGRID ! grid type +CHARACTER(LEN=6), INTENT(IN) :: HREC ! name of the parameter +INTEGER, INTENT(IN) :: KFIELD ! integer scalar for complete grid +INTEGER, INTENT(OUT):: KFIELD_EXTEND ! integer scalar for splitted grid +!* 0.2 Declarations of local variables +! +INTEGER :: IIB, IIE, IJB, IJE +!------------------------------------------------------------------------------- +! +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +! +IF (HREC=='IMAX') KFIELD_EXTEND = IIE-IIB+1 + 2*NHALO +IF (HREC=='JMAX') KFIELD_EXTEND = IJE-IJB+1 + 2*NHALO +! +!------------------------------------------------------------------------------- +END SUBROUTINE EXTEND_GRID_PARAMETERN0_MNH diff --git a/src/MNH/goto_model_surfex_mnh.f90 b/src/MNH/goto_model_surfex_mnh.f90 new file mode 100644 index 000000000..818bd1026 --- /dev/null +++ b/src/MNH/goto_model_surfex_mnh.f90 @@ -0,0 +1,168 @@ +!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +!--------------- special set of characters for RCS information +!----------------------------------------------------------------- +! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/goto_model_surfex_mnh.f90 +!----------------------------------------------------------------- +!####################### +MODULE MODI_GOTO_MODEL_SURFEX_MNH + !####################### + ! + INTERFACE + ! ############################### + SUBROUTINE GOTO_MODEL_SURFEX_MNH(KMI, KINFO_ll) + ! ############################### + !! + !! PURPOSE + !! ------- + !! + !! Initializes local sizes in SURFEX module MODD_SURF_ATM_n for model KMI + !! and calls GOTO_MODEL(KMI) + !! GOTO_MODEL_ll(KMI, KINFO_ll) + !! + !! METHOD + !! ------ + !! + !! EXTERNAL + !! -------- + !! + !! + !! IMPLICIT ARGUMENTS + !! ------------------ + !! + !! + !! REFERENCE + !! --------- + !! + !! AUTHOR + !! ------ + !! + !! M. Moge LA - CNRS + !! + !! MODIFICATION + !! ------------ + !! + !! Original 08/2015 + !---------------------------------------------------------------------------- + ! + !* 0. DECLARATION + ! ----------- + ! +! USE MODE_ll + USE MODE_TOOLS_ll, ONLY : GET_GLOBALDIMS_ll, GET_DIM_PHYS_ll + USE MODD_SURF_ATM_n, ONLY : NDIM_FULL, NSIZE_FULL, NIMAX_SURF_ll, NJMAX_SURF_ll, NIMAX_SURF_LOC, NJMAX_SURF_LOC +! USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT + ! +! USE MODI_GOTO_SURFEX + ! + IMPLICIT NONE + ! + !* 0.1 Declaration of dummy arguments + ! ------------------------------ + ! + INTEGER, INTENT(IN) :: KMI !model id + INTEGER, INTENT(OUT) :: KINFO_ll + END SUBROUTINE GOTO_MODEL_SURFEX_MNH + ! + END INTERFACE + ! +END MODULE MODI_GOTO_MODEL_SURFEX_MNH +! ############################### + SUBROUTINE GOTO_MODEL_SURFEX_MNH(KMI, KINFO_ll) +! ############################### +!! +!! PURPOSE +!! ------- +!! +!! Initializes local sizes in SURFEX module MODD_SURF_ATM_n for model KMI +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M. Moge LA - CNRS +!! +!! MODIFICATION +!! ------------ +!! +!! Original 08/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODD_SURF_ATM_n, ONLY : NDIM_FULL, NSIZE_FULL, NIMAX_SURF_ll, NJMAX_SURF_ll, NIMAX_SURF_LOC, NJMAX_SURF_LOC +! +!USE MODE_ll +USE MODE_TOOLS_ll, ONLY : GET_GLOBALDIMS_ll, GET_DIM_PHYS_ll +USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT +USE MODD_VAR_ll, ONLY : YSPLITTING +USE MODE_MODELN_HANDLER +! +USE MODI_GOTO_SURFEX +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! +INTEGER, INTENT(IN) :: KMI !model id +INTEGER, INTENT(OUT) :: KINFO_ll +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +!INTEGER :: IINFO_ll ! return code of // routines +INTEGER :: IMI ! return code of // routines +CHARACTER*1 :: HSPLIT +! +!------------------------------------------------------------------------------ +! +!IMI = GET_CURRENT_MODEL_INDEX() +!IF ( KMI /= IMI ) THEN +! WRITE(*,*) "ERROR IN PGD_GOTO_MODEL_SURFEX_MNH : current MNH model ", IMI, " is different from target SURFEX model ", KMI +! CALL ABORT +!ENDIF +!! il faudrait faire ce test mais cette GOTO_SURFEX ne marche pas, donc on laisse tomber et on bricole +! +!------------------------------------------------------------------------------ +! +!CALL GOTO_SURFEX(KMI, LKFROM) ! cette routine plante, donc on ne touche pas aux modeles +!CALL GO_TOMODEL_ll(KMI,KINFO_ll) +!CALL GOTO_MODEL(KMI) +! +CALL GET_GLOBALDIMS_ll(NIMAX_SURF_ll,NJMAX_SURF_ll) +NDIM_FULL = NIMAX_SURF_ll*NJMAX_SURF_ll +! +IF ( YSPLITTING == "BSPLITTING" ) THEN + HSPLIT = 'B' +ELSE IF ( YSPLITTING == "XSPLITTING" ) THEN + HSPLIT = 'X' +ELSE IF ( YSPLITTING == "YSPLITTING" ) THEN + HSPLIT = 'Y' +ELSE + HSPLIT = '' +ENDIF +CALL GET_DIM_PHYS_ll(HSPLIT,NIMAX_SURF_LOC,NJMAX_SURF_LOC) +NSIZE_FULL = NIMAX_SURF_LOC*NJMAX_SURF_LOC +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE GOTO_MODEL_SURFEX_MNH diff --git a/src/MNH/mnh_surf_grid_io_init.f90 b/src/MNH/mnh_surf_grid_io_init.f90 new file mode 100644 index 000000000..4d3712227 --- /dev/null +++ b/src/MNH/mnh_surf_grid_io_init.f90 @@ -0,0 +1,158 @@ +!----------------------------------------------------------------- +!--------------- special set of characters for RCS information +!----------------------------------------------------------------- +! $Source$ $Revision$ +! masdev4_7 BUG1 2007/06/15 17:47:18 +!----------------------------------------------------------------- +!####################### +MODULE MODI_MNH_SURF_GRID_IO_INIT + !####################### + ! + INTERFACE + ! ############################### + SUBROUTINE MNH_SURF_GRID_IO_INIT(KIMAX,KJMAX) + ! ############################### + !! + !! PURPOSE + !! ------- + !! + !! Initializes parallel data structures for I/O in surfex (for grid reading) + !! + !! METHOD + !! ------ + !! + !! EXTERNAL + !! -------- + !! + !! + !! IMPLICIT ARGUMENTS + !! ------------------ + !! + !! + !! REFERENCE + !! --------- + !! + !! AUTHOR + !! ------ + !! + !! M.Moge CNRS - LA + !! + !! MODIFICATION + !! ------------ + !! + !! Original 19/03/2015 + !---------------------------------------------------------------------------- + ! + !* 0. DECLARATION + ! ----------- + ! + USE MODE_ll + USE MODE_FM + USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT, JPMODELMAX + USE MODD_CONF, ONLY : CPROGRAM, L1D, L2D, LPACK + ! + USE MODE_SPLITTINGZ_ll + ! + USE MODI_GET_SURF_GRID_DIM_N + USE MODI_GET_LUOUT + ! + IMPLICIT NONE + ! + !* 0.1 Declaration of dummy arguments + ! ------------------------------ + ! + INTEGER, INTENT(IN) :: KIMAX ! number of points in X direction + INTEGER, INTENT(IN) :: KJMAX ! number of points in Y direction + END SUBROUTINE MNH_SURF_GRID_IO_INIT + ! + END INTERFACE +END MODULE MODI_MNH_SURF_GRID_IO_INIT +! ############################### + SUBROUTINE MNH_SURF_GRID_IO_INIT(KIMAX,KJMAX) +! ############################### +!! +!! PURPOSE +!! ------- +!! +!! Initializes parallel data structures for I/O in surfex (for grid reading) +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M.Moge CNRS - LA +!! +!! MODIFICATION +!! ------------ +!! +!! Original 19/03/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODE_ll +USE MODE_FM +USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT, JPMODELMAX +USE MODD_CONF, ONLY : CPROGRAM, L1D, L2D, LPACK +! +!JUANZ +USE MODE_SPLITTINGZ_ll +!JUANZ +! +USE MODI_GET_SURF_GRID_DIM_N +USE MODI_GET_LUOUT +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! +INTEGER, INTENT(IN) :: KIMAX ! number of points in X direction +INTEGER, INTENT(IN) :: KJMAX ! number of points in Y direction +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +INTEGER :: IINFO_ll ! return code of // routines +INTEGER :: ILUOUT ! output listing logical unit +! +!------------------------------------------------------------------------------ +! +IF (CPROGRAM=='IDEAL ' .OR. CPROGRAM=='SPAWN ' .OR. CPROGRAM=='REAL ') RETURN +! +L1D=(KIMAX==1).AND.(KJMAX==1) +L2D=(KIMAX/=1).AND.(KJMAX==1) +LPACK=L1D.OR.L2D +CALL SET_FMPACK_ll(L1D,L2D,LPACK) +CALL SET_JP_ll(JPMODELMAX,JPHEXT,JPVEXT,JPHEXT) +CALL SET_DAD0_ll() +CALL SET_DIM_ll(KIMAX, KJMAX, 1) +CALL SET_LBX_ll('OPEN',1) +CALL SET_LBY_ll('OPEN', 1) +CALL SET_XRATIO_ll(1, 1) +CALL SET_YRATIO_ll(1, 1) +CALL SET_XOR_ll(1, 1) +CALL SET_XEND_ll(KIMAX+2*JPHEXT, 1) +CALL SET_YOR_ll(1, 1) +CALL SET_YEND_ll(KJMAX+2*JPHEXT, 1) +CALL SET_DAD_ll(0, 1) +CALL INI_PARAZ_ll(IINFO_ll) +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE MNH_SURF_GRID_IO_INIT diff --git a/src/SURFEX/extend_grid_on_halo.F90 b/src/SURFEX/extend_grid_on_halo.F90 new file mode 100644 index 000000000..a7d7a0550 --- /dev/null +++ b/src/SURFEX/extend_grid_on_halo.F90 @@ -0,0 +1,87 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +!! ########################################################### + SUBROUTINE EXTEND_GRID_ON_HALO(HPROGRAM,KGRID_PAR,PGRID_PAR) +! ########################################################### +!! +!! PURPOSE +!! ------- +!! This program extends a splitted PGD grid on the SURFEX halo +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M.Moge CNRS - LA +!! +!! MODIFICATION +!! ------------ +!! +!! Original 01/03/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODD_SURF_ATM_GRID_n, ONLY : CGRID!, XGRID_PAR, NGRID_PAR +USE MODD_SURF_ATM_n, ONLY : NDIM_FULL, NSIZE_FULL +! +USE MODI_EXTEND_GRID_ON_HALO_CONF_PROJ +USE MODI_EXTEND_GRID_ON_HALO_CARTESIAN +USE MODI_GET_SIZE_FULL_n +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling +INTEGER, INTENT(INOUT) :: KGRID_PAR ! size of PGRID_PAR pointer +REAL, DIMENSION(:), POINTER, INTENT(INOUT) :: PGRID_PAR ! parameters defining this grid +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE + CHARACTER(LEN=100) :: YCOMMENT +INTEGER :: IRESP ! error return code +INTEGER :: IHALO +!------------------------------------------------------------------------------ +IF (LHOOK) CALL DR_HOOK('SPLIT_GRID',0,ZHOOK_HANDLE) +! +SELECT CASE(CGRID) + + CASE('CONF PROJ ') + CALL EXTEND_GRID_ON_HALO_CONF_PROJ(HPROGRAM,NDIM_FULL,NSIZE_FULL,KGRID_PAR,PGRID_PAR) + CASE('CARTESIAN ') + CALL EXTEND_GRID_ON_HALO_CARTESIAN(HPROGRAM,NDIM_FULL,NSIZE_FULL,KGRID_PAR,PGRID_PAR) + CASE DEFAULT + CALL GET_SIZE_FULL_n(HPROGRAM,NDIM_FULL,NSIZE_FULL) + +END SELECT +! + +IF (LHOOK) CALL DR_HOOK('SPLIT_GRID',1,ZHOOK_HANDLE) +!_______________________________________________________________________________ +! +END SUBROUTINE EXTEND_GRID_ON_HALO diff --git a/src/SURFEX/extend_grid_on_halo_cartesian.F90 b/src/SURFEX/extend_grid_on_halo_cartesian.F90 new file mode 100644 index 000000000..203e5181d --- /dev/null +++ b/src/SURFEX/extend_grid_on_halo_cartesian.F90 @@ -0,0 +1,117 @@ + +! ########################################################### + SUBROUTINE EXTEND_GRID_ON_HALO_CARTESIAN(HPROGRAM,KDIM_FULL,KSIZE_FULL,KGRID_PAR,PGRID_PAR) +! ########################################################### +!! +!! PURPOSE +!! ------- +!! This program extends a splitted PGD grid on the SURFEX halo +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M.Moge CNRS - LA +!! +!! MODIFICATION +!! ------------ +!! +!! Original 01/03/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODE_GRIDTYPE_CARTESIAN +USE MODE_EXTEND_GRID_PARAMETER +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! host program +INTEGER, INTENT(IN) :: KDIM_FULL ! total number of points +INTEGER, INTENT(OUT) :: KSIZE_FULL! number of points on this processor +INTEGER, INTENT(INOUT) :: KGRID_PAR ! size of PGRID_PAR pointer +REAL, DIMENSION(:), POINTER :: PGRID_PAR ! parameters defining this grid +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!* original grid +REAL :: ZLAT0, ZLON0 +INTEGER :: IIMAX, IJMAX +REAL, DIMENSION(KDIM_FULL) :: ZX, ZY, ZDX, ZDY +! +!* extended grid +INTEGER :: IIMAX_EXTENDED, IJMAX_EXTENDED +REAL, DIMENSION(:), ALLOCATABLE :: ZX_EXTENDED, ZY_EXTENDED, ZDX_EXTENDED, ZDY_EXTENDED +! +!------------------------------------------------------------------------------ +IF (LHOOK) CALL DR_HOOK('EXTEND_GRID_CARTESIAN',0,ZHOOK_HANDLE) +! +!* 1. Gets Parameters of the Grid +! + CALL GET_GRIDTYPE_CARTESIAN(PGRID_PAR,ZLAT0,ZLON0, & + IIMAX,IJMAX, & + ZX,ZY,ZDX,ZDY ) +! +! +!* 2. Splits the (pertinent) parameters of the grid +! + CALL EXTEND_GRID_PARAMETERN0(HPROGRAM,'CARTESIAN ','IMAX ',IIMAX,IIMAX_EXTENDED) + CALL EXTEND_GRID_PARAMETERN0(HPROGRAM,'CARTESIAN ','JMAX ',IJMAX,IJMAX_EXTENDED) +! +KSIZE_FULL = IIMAX_EXTENDED * IJMAX_EXTENDED +! +ALLOCATE(ZX_EXTENDED (KSIZE_FULL)) +ALLOCATE(ZY_EXTENDED (KSIZE_FULL)) +ALLOCATE(ZDX_EXTENDED(KSIZE_FULL)) +ALLOCATE(ZDY_EXTENDED(KSIZE_FULL)) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CARTESIAN ','XX ',SIZE(ZX),KSIZE_FULL,IIMAX,IJMAX,ZX,ZX_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CARTESIAN ','YY ',SIZE(ZY),KSIZE_FULL,IIMAX,IJMAX,ZY,ZY_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CARTESIAN ','DX ',SIZE(ZDX),KSIZE_FULL,IIMAX,IJMAX,ZDX,ZDX_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CARTESIAN ','DY ',SIZE(ZDY),KSIZE_FULL,IIMAX,IJMAX,ZDY,ZDY_EXTENDED) +! +! +!* 3. Stores Parameters of the Grid in grid pointer +! +NULLIFY(PGRID_PAR) + CALL PUT_GRIDTYPE_CARTESIAN(PGRID_PAR,ZLAT0,ZLON0, & + IIMAX_EXTENDED,IJMAX_EXTENDED, & + ZX_EXTENDED,ZY_EXTENDED,ZDX_EXTENDED,ZDY_EXTENDED ) + ! +! +KGRID_PAR = SIZE(PGRID_PAR) +! +DEALLOCATE(ZX_EXTENDED ) +DEALLOCATE(ZY_EXTENDED ) +DEALLOCATE(ZDX_EXTENDED) +DEALLOCATE(ZDY_EXTENDED) +! +IF (LHOOK) CALL DR_HOOK('EXTEND_GRID_CARTESIAN',1,ZHOOK_HANDLE) +!_______________________________________________________________________________ +! +END SUBROUTINE EXTEND_GRID_ON_HALO_CARTESIAN diff --git a/src/SURFEX/extend_grid_on_halo_conf_proj.F90 b/src/SURFEX/extend_grid_on_halo_conf_proj.F90 new file mode 100644 index 000000000..517221e6a --- /dev/null +++ b/src/SURFEX/extend_grid_on_halo_conf_proj.F90 @@ -0,0 +1,117 @@ + +! ########################################################### + SUBROUTINE EXTEND_GRID_ON_HALO_CONF_PROJ(HPROGRAM,KDIM_FULL,KSIZE_FULL,KGRID_PAR,PGRID_PAR) +! ########################################################### +!! +!! PURPOSE +!! ------- +!! This program extends a splitted PGD grid on the SURFEX halo +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M.Moge CNRS - LA +!! +!! MODIFICATION +!! ------------ +!! +!! Original 01/03/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODE_GRIDTYPE_CONF_PROJ +USE MODE_EXTEND_GRID_PARAMETER +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! host program +INTEGER, INTENT(IN) :: KDIM_FULL ! total number of points +INTEGER, INTENT(OUT) :: KSIZE_FULL! number of points on this processor +INTEGER, INTENT(INOUT) :: KGRID_PAR ! size of PGRID_PAR pointer +REAL, DIMENSION(:), POINTER :: PGRID_PAR ! parameters defining this grid +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!* original grid +REAL :: ZLAT0, ZLON0, ZRPK, ZBETA, ZLATOR, ZLONOR +INTEGER :: IIMAX, IJMAX +REAL, DIMENSION(PGRID_PAR(11)) :: ZX, ZY, ZDX, ZDY +! +!* extended grid +INTEGER :: IIMAX_EXTENDED, IJMAX_EXTENDED +REAL, DIMENSION(:), ALLOCATABLE :: ZX_EXTENDED, ZY_EXTENDED, ZDX_EXTENDED, ZDY_EXTENDED +! +!------------------------------------------------------------------------------ +IF (LHOOK) CALL DR_HOOK('EXTEND_GRID_CONF_PROJ',0,ZHOOK_HANDLE) +! +!* 1. Gets Parameters of the Grid +! + CALL GET_GRIDTYPE_CONF_PROJ(PGRID_PAR,ZLAT0,ZLON0,ZRPK,ZBETA,& + ZLATOR,ZLONOR,IIMAX,IJMAX, & + ZX,ZY,ZDX,ZDY ) +! +! +!* 2. Splits the (pertinent) parameters of the grid +! + CALL EXTEND_GRID_PARAMETERN0(HPROGRAM,'CONF PROJ ','IMAX ',IIMAX,IIMAX_EXTENDED) + CALL EXTEND_GRID_PARAMETERN0(HPROGRAM,'CONF PROJ ','JMAX ',IJMAX,IJMAX_EXTENDED) +! +KSIZE_FULL = IIMAX_EXTENDED * IJMAX_EXTENDED +! +ALLOCATE(ZX_EXTENDED (KSIZE_FULL)) +ALLOCATE(ZY_EXTENDED (KSIZE_FULL)) +ALLOCATE(ZDX_EXTENDED(KSIZE_FULL)) +ALLOCATE(ZDY_EXTENDED(KSIZE_FULL)) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CONF PROJ ','XX ',SIZE(ZX),KSIZE_FULL,IIMAX,IJMAX,ZX,ZX_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CONF PROJ ','YY ',SIZE(ZY),KSIZE_FULL,IIMAX,IJMAX,ZY,ZY_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CONF PROJ ','DX ',SIZE(ZDX),KSIZE_FULL,IIMAX,IJMAX,ZDX,ZDX_EXTENDED) + CALL EXTEND_GRID_PARAMETERX1(HPROGRAM,'CONF PROJ ','DY ',SIZE(ZDY),KSIZE_FULL,IIMAX,IJMAX,ZDY,ZDY_EXTENDED) +! +! +!* 3. Stores Parameters of the Grid in grid pointer +! +NULLIFY(PGRID_PAR) + CALL PUT_GRIDTYPE_CONF_PROJ(PGRID_PAR,ZLAT0,ZLON0,ZRPK,ZBETA, & + ZLATOR,ZLONOR,IIMAX_EXTENDED,IJMAX_EXTENDED, & + ZX_EXTENDED,ZY_EXTENDED,ZDX_EXTENDED,ZDY_EXTENDED ) + ! +! +KGRID_PAR = SIZE(PGRID_PAR) +! +DEALLOCATE(ZX_EXTENDED ) +DEALLOCATE(ZY_EXTENDED ) +DEALLOCATE(ZDX_EXTENDED) +DEALLOCATE(ZDY_EXTENDED) +! +IF (LHOOK) CALL DR_HOOK('EXTEND_GRID_CONF_PROJ',1,ZHOOK_HANDLE) +!_______________________________________________________________________________ +! +END SUBROUTINE EXTEND_GRID_ON_HALO_CONF_PROJ diff --git a/src/SURFEX/goto_model_mnh.F90 b/src/SURFEX/goto_model_mnh.F90 new file mode 100644 index 000000000..6379e26d8 --- /dev/null +++ b/src/SURFEX/goto_model_mnh.F90 @@ -0,0 +1,141 @@ +!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +!--------------- special set of characters for RCS information +!----------------------------------------------------------------- +! $Source: /home/cvsroot/MNH-VX-Y-Z/src/SURFEX/goto_model_surfex_mnh.F90 +!----------------------------------------------------------------- +!####################### +MODULE MODI_GOTO_MODEL_MNH + !####################### + ! + INTERFACE + ! ############################### + SUBROUTINE GOTO_MODEL_MNH(HPROGRAM, KMI, KINFO_ll) + ! ############################### + !! + !! PURPOSE + !! ------- + !! + !! Calls GOTO_MODEL_SURFEX_MNH to + !! initialize local sizes in SURFEX module MODD_SURF_ATM_n for model KMI + !! and call GOTO_MODEL(KMI) + !! GOTO_MODEL_ll(KMI, KINFO_ll) + !! + !! METHOD + !! ------ + !! + !! EXTERNAL + !! -------- + !! + !! + !! IMPLICIT ARGUMENTS + !! ------------------ + !! + !! + !! REFERENCE + !! --------- + !! + !! AUTHOR + !! ------ + !! + !! M. Moge LA - CNRS + !! + !! MODIFICATION + !! ------------ + !! + !! Original 08/2015 + !---------------------------------------------------------------------------- + ! + !* 0. DECLARATION + ! ----------- + ! + IMPLICIT NONE + ! + !* 0.1 Declaration of dummy arguments + ! ------------------------------ + ! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program + INTEGER, INTENT(IN) :: KMI !model id + INTEGER, INTENT(OUT) :: KINFO_ll + END SUBROUTINE GOTO_MODEL_MNH + ! + END INTERFACE + ! +END MODULE MODI_GOTO_MODEL_MNH +! ############################### + SUBROUTINE GOTO_MODEL_MNH(HPROGRAM, KMI, KINFO_ll) +! ############################### +!! +!! PURPOSE +!! ------- +!! +!! Calls GOTO_MODEL_SURFEX_MNH to +!! initialize local sizes in SURFEX module MODD_SURF_ATM_n for model KMI +!! and call GOTO_MODEL(KMI) +!! GOTO_MODEL_ll(KMI, KINFO_ll) +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M. Moge LA - CNRS +!! +!! MODIFICATION +!! ------------ +!! +!! Original 08/2015 +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +! +#ifdef MNH +USE MODI_GOTO_MODEL_SURFEX_MNH +#endif +! +IMPLICIT NONE +! +!* 0.1 Declaration of dummy arguments +! ------------------------------ +! +CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program +INTEGER, INTENT(IN) :: KMI !model id +INTEGER, INTENT(OUT) :: KINFO_ll +! +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +INTEGER :: IMI ! return code of // routines +CHARACTER*1 :: HSPLIT +! +!------------------------------------------------------------------------------ +! +IF (HPROGRAM=='MESONH') THEN +#ifdef MNH + CALL GOTO_MODEL_SURFEX_MNH(KMI, KINFO_ll) +#else + KINFO_ll = 0 +#endif +ENDIF +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE GOTO_MODEL_MNH diff --git a/src/SURFEX/hor_interpol_1cov.F90 b/src/SURFEX/hor_interpol_1cov.F90 new file mode 100644 index 000000000..8273f3bdc --- /dev/null +++ b/src/SURFEX/hor_interpol_1cov.F90 @@ -0,0 +1,133 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL * - Call the interpolation of a surface field +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +USE MODD_PREP, ONLY : CINGRID_TYPE, CINTERP_TYPE +! +USE MODI_HOR_INTERPOL_GAUSS_1COV +USE MODI_HOR_INTERPOL_ROTLATLON_1COV +USE MODI_HOR_INTERPOL_AROME_1COV +USE MODI_HOR_INTERPOL_CONF_PROJ_1COV +USE MODI_HOR_INTERPOL_CARTESIAN_1COV +USE MODI_HOR_INTERPOL_LATLON_1COV +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +USE MODI_ABOR1_SFX +! +USE MODI_HOR_INTERPOL_BUFFER +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +INTEGER :: JL ! loop counter +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------------- +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_1COV',0,ZHOOK_HANDLE) +SELECT CASE (CINTERP_TYPE) +! +!* 1. Interpolation with horibl (from gaussian, Legendre or regular grid) +! ------------------------------------------------------------------- +! + CASE('HORIBL') + SELECT CASE(CINGRID_TYPE) +! +!* 1.1 Interpolation from gaussian or Legendre +! + CASE ('GAUSS ') + CALL HOR_INTERPOL_GAUSS_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! +!* 1.2 Interpolation from regular grid +! + CASE ('AROME ') + CALL HOR_INTERPOL_AROME_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! +!* 1.3 Interpolation from regular lat/lon coord +! + CASE ('LATLON ') + CALL HOR_INTERPOL_LATLON_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! +!* 1.4 Interpolation from rotated lat/lon coord +! + CASE ('ROTLATLON ') + CALL HOR_INTERPOL_ROTLATLON_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! + CASE DEFAULT + CALL ABOR1_SFX('HOR_INTERPOL_1COV: WRONG GRID TYPE'//CINGRID_TYPE) + + END SELECT +! +!* 2. Prescribed uniform field +! ------------------------ +! + CASE('UNIF ') + PFIELDOUT(:) = PFIELDIN(1) +! +!* 3. Bilinear interpolation +! ---------------------- +! + CASE('BILIN ') + SELECT CASE(CINGRID_TYPE) + CASE ('CONF PROJ ') + CALL HOR_INTERPOL_CONF_PROJ_1COV(KLUOUT,PFIELDIN,PFIELDOUT) + CASE ('CARTESIAN ') + CALL HOR_INTERPOL_CARTESIAN_1COV(KLUOUT,PFIELDIN,PFIELDOUT) + END SELECT +! +!* 4. no interpolation, only packing +! ------------------------------ +! + CASE('BUFFER') + CALL HOR_INTERPOL_BUFFER_1COV(KLUOUT,PFIELDIN,PFIELDOUT) + +! +!* 4. no interpolation +! ---------------- +! + CASE('NONE ') + PFIELDOUT(:) = PFIELDIN(:) + + CASE DEFAULT + CALL ABOR1_SFX('HOR_INTERPOL_1COV: WRONG INTERPOLATION TYPE'//CINTERP_TYPE) + +END SELECT +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_1COV',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_1COV diff --git a/src/SURFEX/hor_interpol_arome_1cov.F90 b/src/SURFEX/hor_interpol_arome_1cov.F90 new file mode 100644 index 000000000..3ea0f1cfe --- /dev/null +++ b/src/SURFEX/hor_interpol_arome_1cov.F90 @@ -0,0 +1,105 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_AROME_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_AROME * - Interpolation from an AROME grid on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_AROME) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT,LINTERP +USE MODD_GRID_AROME, ONLY : XX, XY, NX, NY, XLAT0, XLON0, XLATOR, XLONOR, XRPK, XBETA +USE MODD_GRID_GRIB, ONLY : NNI +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODE_GRIDTYPE_CONF_PROJ +USE MODI_HORIBL_SURF +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +REAL, DIMENSION(:), ALLOCATABLE :: ZX ! X coordinate +REAL, DIMENSION(:), ALLOCATABLE :: ZY ! Y coordinate +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKIN ! input mask +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKOUT ! output mask +INTEGER :: INO ! output number of points +INTEGER :: JL ! loop counter +INTEGER, DIMENSION(:), ALLOCATABLE :: IX ! number of points on each line +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------------- +! +!* 1. Allocations +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_AROME_1COV',0,ZHOOK_HANDLE) +INO = SIZE(XLAT_OUT) +! +ALLOCATE(IMASKIN (NNI)) +! +ALLOCATE(ZX (INO)) +ALLOCATE(ZY (INO)) +ALLOCATE(IMASKOUT(INO)) +IMASKOUT = 1 +ALLOCATE(IX(NY)) +IX=NX +! +!* 2. Transformation of latitudes/longitudes into metric coordinates of input grid +! + CALL XY_CONF_PROJ(XLAT0,XLON0,XRPK,XBETA,XLATOR,XLONOR,ZX,ZY,XLAT_OUT,XLON_OUT) +! +! +!* 3. Input mask +! +IMASKIN(:) = 1 +WHERE(PFIELDIN(:)==XUNDEF) IMASKIN = 0 +! +! +!* 4. Interpolation with horibl +! +CALL HORIBL_SURF(0.,0.,XY,XX,NY,IX,NNI,PFIELDIN(:),INO,ZX,ZY,PFIELDOUT(:), & + .FALSE.,KLUOUT,LINTERP,IMASKIN,IMASKOUT) +! +!* 5. Deallocations +! +DEALLOCATE(ZX) +DEALLOCATE(IX) +DEALLOCATE(ZY) +DEALLOCATE(IMASKIN ) +DEALLOCATE(IMASKOUT) +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_AROME_1COV',1,ZHOOK_HANDLE) + +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_AROME_1COV diff --git a/src/SURFEX/hor_interpol_buffer_1cov.F90 b/src/SURFEX/hor_interpol_buffer_1cov.F90 new file mode 100644 index 000000000..e371eb79e --- /dev/null +++ b/src/SURFEX/hor_interpol_buffer_1cov.F90 @@ -0,0 +1,78 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_BUFFER_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_BUFFER_1COV * - Only extrapolation on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_BUFFER) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : CMASK +! +USE MODD_GRID_BUFFER, ONLY : NNI +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODI_PACK_SAME_RANK +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +USE MODI_GET_SURF_MASK_n +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKOUT ! output mask +INTEGER :: INO ! output number of points +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------ +! +!* 1. Initialisation of the output mask +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_BUFFER_1COV',0,ZHOOK_HANDLE) +INO = SIZE(PFIELDOUT,1) +ALLOCATE(IMASKOUT(INO)) + CALL GET_SURF_MASK_n(CMASK,INO,IMASKOUT,NNI,KLUOUT) +! +!* 2. Mask the input field with the output mask +!!mask du tableau de taille FULL en fonction du type de surface + CALL PACK_SAME_RANK(IMASKOUT,PFIELDIN,PFIELDOUT) +! +!* 6. Deallocations +! +DEALLOCATE(IMASKOUT) +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_BUFFER_1COV',1,ZHOOK_HANDLE) + +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_BUFFER_1COV diff --git a/src/SURFEX/hor_interpol_cartesian_1cov.F90 b/src/SURFEX/hor_interpol_cartesian_1cov.F90 new file mode 100644 index 000000000..141a06205 --- /dev/null +++ b/src/SURFEX/hor_interpol_cartesian_1cov.F90 @@ -0,0 +1,107 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_CARTESIAN_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_CARTESIAN_1COV * - Interpolation from a CARTESIAN grid on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_CARTESIAN) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! ################################################################################# +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT, XX_OUT, XY_OUT, LINTERP +USE MODD_GRID_CARTESIAN, ONLY : XX, XY, NX, NY +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODE_GRIDTYPE_CARTESIAN +USE MODI_BILIN +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +REAL, DIMENSION(:), ALLOCATABLE :: ZX ! X coordinate +REAL, DIMENSION(:), ALLOCATABLE :: ZY ! Y coordinate +INTEGER :: INO ! output number of points +REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELDIN ! input field +! +INTEGER :: JI ! loop index +INTEGER :: JJ ! loop index +INTEGER :: JL ! loop index +REAL(KIND=JPRB) :: ZHOOK_HANDLE + +!------------------------------------------------------------------------------------- +! +!* 1. Allocations +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_CARTESIAN_1COV',0,ZHOOK_HANDLE) +INO = SIZE(XX_OUT) +! +ALLOCATE(ZX (INO)) +ALLOCATE(ZY (INO)) +! +!* 2. Transformation of latitudes/longitudes into metric coordinates of output grid +! +!* WARNING : here, because the input grid is not geographic, one assumes that +! coordinates are coherent between input and output grid, but without +! any way to check it. +! +ZX = XX_OUT +ZY = XY_OUT +! +! +!* 3. Put input field on its squared grid +! +ALLOCATE(ZFIELDIN(NX,NY)) +! +DO JJ=1,NY + DO JI=1,NX + ZFIELDIN(JI,JJ) = PFIELDIN(JI+NX*(JJ-1)) + END DO +END DO +! +!* 4. Interpolation with bilinear +! +CALL BILIN(KLUOUT,XX,XY,ZFIELDIN(:,:),ZX,ZY,PFIELDOUT(:),LINTERP) +! +! +!* 5. Deallocations +! +! +DEALLOCATE(ZX) +DEALLOCATE(ZY) +DEALLOCATE(ZFIELDIN) +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_CARTESIAN_1COV',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_CARTESIAN_1COV diff --git a/src/SURFEX/hor_interpol_conf_proj_1cov.F90 b/src/SURFEX/hor_interpol_conf_proj_1cov.F90 new file mode 100644 index 000000000..fc75d8ac0 --- /dev/null +++ b/src/SURFEX/hor_interpol_conf_proj_1cov.F90 @@ -0,0 +1,142 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_CONF_PROJ_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +!! +!! PURPOSE +!! ------- same as HOR_INTERPOL_CONF_PROJ on only 1 vertical level +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! MODIFICATION +!! ------------ +!! +!! 02/04/12 M. Tomasini Add an index in the second dimension of the ISBA +!! array variables for BILIN interpolation routine to +!! not bug in case 2D (this is not the more beautiful +!! method; the BILIN routine should better be adapted) +!! Search ! Ajout MT +!! 10/02/15 M.Moge using SIZE(PFIELDOUT,1) instead of SIZE(XLAT_OUT) +!------------------------------------------------------------------------------- +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT, LINTERP +USE MODD_GRID_CONF_PROJ, ONLY : XX, XY, NX, NY, XLAT0, XLON0, XLATORI, & + XLONORI, XRPK, XBETA +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODE_GRIDTYPE_CONF_PROJ +USE MODI_BILIN +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +REAL, DIMENSION(:), ALLOCATABLE :: ZX,ZY ! coordinate of the output field +REAL, DIMENSION(:), ALLOCATABLE :: ZX_DUPLIQUE ! X coordinate of the output field ! Ajout MT +REAL, DIMENSION(:), ALLOCATABLE :: ZY_DUPLIQUE ! Y coordinate of the output field ! Ajout MT +REAL, DIMENSION(:), ALLOCATABLE :: ZXY_DUPLIQUE ! Y coordinate of the input field ! Ajout MT +REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELDIN ! input field +REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELDIN_DUPLIQUE ! input field ! Ajout MT +REAL, DIMENSION(:), ALLOCATABLE :: ZFIELDOUT_DUPLIQUE ! interpolated output field ! Ajout MT +! +INTEGER :: INO ! output number of points +INTEGER :: JI,JJ,JL ! loop index +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +LOGICAL, DIMENSION(:), ALLOCATABLE :: GINTERP_DUPLIQUE ! .true. where physical value is needed ! Ajout MT +!------------------------------------------------------------------------------------- +! +!* 1. Allocations +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_CONF_PROJ_1COV',0,ZHOOK_HANDLE) +INO = SIZE(PFIELDOUT) +! +ALLOCATE(ZX (INO)) +ALLOCATE(ZY (INO)) +! +IF (NY==1) THEN ! Ajout MT + ALLOCATE(ZXY_DUPLIQUE(2),ZFIELDIN_DUPLIQUE(NX,2)) + ALLOCATE(ZX_DUPLIQUE(2*INO),ZY_DUPLIQUE(2*INO),ZFIELDOUT_DUPLIQUE(2*INO)) + ALLOCATE(GINTERP_DUPLIQUE(SIZE(ZFIELDOUT_DUPLIQUE))) +END IF +! +!* 2. Transformation of latitudes/longitudes into metric coordinates of output grid +! + CALL XY_CONF_PROJ(XLAT0,XLON0,XRPK,XBETA,XLATORI,XLONORI, & + ZX,ZY,XLAT_OUT,XLON_OUT ) +! +!* 3. Put input field on its squared grid +! +ALLOCATE(ZFIELDIN(NX,NY)) +! +DO JJ=1,NY + DO JI=1,NX + ZFIELDIN(JI,JJ) = PFIELDIN(JI+NX*(JJ-1)) + END DO +END DO +! +IF (NY==1) THEN ! Ajout MT + ZFIELDIN_DUPLIQUE(:,1)=ZFIELDIN(:,1) + ZFIELDIN_DUPLIQUE(:,2)=ZFIELDIN(:,1) + ZXY_DUPLIQUE(1)=XY(1) + ZXY_DUPLIQUE(2)=XY(1)+10000. + ZX_DUPLIQUE(1:INO) =ZX(:) + ZX_DUPLIQUE(INO+1:2*INO)=ZX(:) + ZY_DUPLIQUE(1:INO) =ZY(:) + ZY_DUPLIQUE(INO+1:2*INO)=ZY(:)+10000. + GINTERP_DUPLIQUE(1:INO) =LINTERP(1:INO) + GINTERP_DUPLIQUE(INO+1:2*INO)=LINTERP(1:INO) +END IF +! +!* 4. Interpolation with bilinear +! +IF (NY==1) THEN ! Ajout MT + CALL BILIN(KLUOUT,XX,ZXY_DUPLIQUE,ZFIELDIN_DUPLIQUE(:,:), & + ZX_DUPLIQUE,ZY_DUPLIQUE,ZFIELDOUT_DUPLIQUE(:),GINTERP_DUPLIQUE) + + PFIELDOUT(1:INO)=ZFIELDOUT_DUPLIQUE(1:INO) +ELSE + CALL BILIN(KLUOUT,XX,XY,ZFIELDIN(:,:),ZX,ZY,PFIELDOUT(:),LINTERP) +END IF +! +! +!* 5. Deallocations +! +! +DEALLOCATE(ZX,ZY) +DEALLOCATE(ZFIELDIN) +IF (NY==1) DEALLOCATE(ZXY_DUPLIQUE,ZX_DUPLIQUE,ZY_DUPLIQUE, & + ZFIELDIN_DUPLIQUE,ZFIELDOUT_DUPLIQUE,GINTERP_DUPLIQUE) ! Ajout MT +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_CONF_PROJ_1COV',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_CONF_PROJ_1COV diff --git a/src/SURFEX/hor_interpol_gauss_1cov.F90 b/src/SURFEX/hor_interpol_gauss_1cov.F90 new file mode 100644 index 000000000..5baaf77e2 --- /dev/null +++ b/src/SURFEX/hor_interpol_gauss_1cov.F90 @@ -0,0 +1,254 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_GAUSS_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_GAUSS_1COV* - Interpolation from a gaussian grid on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_GAUSS) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT, LINTERP +USE MODD_GRID_GAUSS, ONLY : XILA1, XILO1, XILA2, XILO2, NINLA, NINLO, NILEN, LROTPOLE, & + XLAP, XLOP, XCOEF +USE MODD_GRID_GRIB, ONLY : NNI +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODI_HORIBL_SURF +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +REAL, DIMENSION(:), ALLOCATABLE :: ZLAT ! latitudes +REAL, DIMENSION(:), ALLOCATABLE :: ZLON ! longitudes +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKIN ! input mask +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKOUT ! output mask +INTEGER :: INO ! output number of points +INTEGER :: JL ! loop counter +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------------- +! +!* 1. Allocations +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_GAUSS_1COV',0,ZHOOK_HANDLE) +INO = SIZE(XLAT_OUT) +! +ALLOCATE(IMASKIN (NNI)) +! +ALLOCATE(ZLAT (INO)) +ALLOCATE(ZLON (INO)) +ALLOCATE(IMASKOUT(INO)) +IMASKOUT = 1 +! +!* 2. Is pole rotated? +! +IF (LROTPOLE) THEN +!* transformation of output latitudes, longitudes into rotated coordinates + CALL ARPEGE_STRETCH_A(INO,XLAP,XLOP,XCOEF,XLAT_OUT,XLON_OUT,ZLAT,ZLON) +ELSE + ZLAT = XLAT_OUT + ZLON = XLON_OUT +END IF +! +!* 3. Input mask +! +! +IMASKIN(:) = 1. +WHERE(PFIELDIN(:)==XUNDEF) IMASKIN = 0. +! +! +!* 4. Interpolation with horibl +! +CALL HORIBL_SURF(XILA1,XILO1,XILA2,XILO2,NINLA,NINLO,NNI,PFIELDIN(:),INO,ZLON,ZLAT,PFIELDOUT(:), & + .FALSE.,KLUOUT,LINTERP,IMASKIN,IMASKOUT) +! +! +!* 5. Deallocations +! +DEALLOCATE(ZLAT) +DEALLOCATE(ZLON) +DEALLOCATE(IMASKIN ) +DEALLOCATE(IMASKOUT) +! +!------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_GAUSS_1COV',1,ZHOOK_HANDLE) +CONTAINS +! +! ########################################################################## +SUBROUTINE ARPEGE_STRETCH_A(KN,PLAP,PLOP,PCOEF,PLAR,PLOR,PLAC,PLOC) +! ########################################################################## +!!**** *ARPEGE_STRETCH_A* - Projection to Arpege stretched grid +!! +!! PURPOSE +!! ------- +!! +!! Projection from standard Lat,Lon grid to Arpege stretched grid +!! +!! METHOD +!! ------ +!! +!! The projection is defined in two steps : +!! 1. A rotation to place the stretching pole at the north pole +!! 2. The stretching +!! This routine is a basic implementation of the informations founded in +!! 'Note de travail Arpege n#3' +!! 'Transformation de coordonnees' +!! J.F.Geleyn 1988 +!! This document describes a slightly different transformation in 3 steps. Only the +!! two first steps are to be taken in account (at the time of writing this paper has +!! not been updated). +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! This routine is based on : +!! 'Note de travail ARPEGE' number 3 +!! by J.F. GELEYN (may 1988) +!! +!! AUTHOR +!! ------ +!! +!! V.Bousquet +!! +!! MODIFICATIONS +!! ------------- +!! +!! Original 07/01/1999 +!! V. Masson 01/2004 Externalization of surface +!! +! +! 0. DECLARATIONS +! --------------- +! + USE MODD_CSTS,ONLY : XPI +! + IMPLICIT NONE +! +! 0.1. Declaration of arguments +! ----------------------------- + + INTEGER, INTENT(IN) :: KN ! Number of points to convert + REAL, INTENT(IN) :: PLAP ! Latitude of stretching pole + REAL, INTENT(IN) :: PLOP ! Longitude of stretching pole + REAL, INTENT(IN) :: PCOEF ! Stretching coefficient + REAL, DIMENSION(KN), INTENT(IN) :: PLAR ! Lat. of points + REAL, DIMENSION(KN), INTENT(IN) :: PLOR ! Lon. of points + REAL, DIMENSION(KN), INTENT(OUT) :: PLAC ! Computed pseudo-lat. of points + REAL, DIMENSION(KN), INTENT(OUT) :: PLOC ! Computed pseudo-lon. of points +! + REAL :: ZSINSTRETCHLA ! Sine of stretching point lat. + REAL :: ZSINSTRETCHLO ! Sine of stretching point lon. + REAL :: ZCOSSTRETCHLA ! Cosine of stretching point lat. + REAL :: ZCOSSTRETCHLO ! Cosine of stretching point lon. + REAL :: ZSINLA ! Sine of computed point latitude + REAL :: ZSINLO ! Sine of computed point longitude + REAL :: ZCOSLA ! Cosine of computed point latitude + REAL :: ZCOSLO ! Cosine of computed point longitude + REAL :: ZSINLAS ! Sine of point's pseudo-latitude + REAL :: ZSINLOS ! Sine of point's pseudo-longitude + REAL :: ZCOSLOS ! Cosine of point's pseudo-lon. + REAL :: ZA,ZB,ZD ! Dummy variables used for + REAL :: ZX,ZY ! computations +! + INTEGER :: JLOOP1 ! Dummy loop counter + REAL(KIND=JPRB) :: ZHOOK_HANDLE +! + IF (LHOOK) CALL DR_HOOK('ARPEGE_STRETCH_A',0,ZHOOK_HANDLE) + ZSINSTRETCHLA = SIN(PLAP*XPI/180.) + ZCOSSTRETCHLA = COS(PLAP*XPI/180.) + ZSINSTRETCHLO = SIN(PLOP*XPI/180.) + ZCOSSTRETCHLO = COS(PLOP*XPI/180.) + ! L = longitude (0 = Greenwich, + toward east) + ! l = latitude (90 = N.P., -90 = S.P.) + ! p stands for stretching pole + PLAC(:) = PLAR(:) * XPI / 180. + PLOC(:) = PLOR(:) * XPI / 180. + ! A = 1 + c.c + ZA = 1. + PCOEF*PCOEF + ! B = 1 - c.c + ZB = 1. - PCOEF*PCOEF + DO JLOOP1=1, KN + ZSINLA = SIN(PLAC(JLOOP1)) + ZCOSLA = COS(PLAC(JLOOP1)) + ZSINLO = SIN(PLOC(JLOOP1)) + ZCOSLO = COS(PLOC(JLOOP1)) + ! X = cos(Lp-L) + ZX = ZCOSLO*ZCOSSTRETCHLO + ZSINLO*ZSINSTRETCHLO + ! Y = sin(Lp-L) + ZY = ZSINSTRETCHLO*ZCOSLO - ZSINLO*ZCOSSTRETCHLO + ! D = (1+c.c) + (1-c.c)(sin lp.sin l + cos lp.cos l.cos(Lp-L)) + ZD = ZA + ZB*(ZSINSTRETCHLA*ZSINLA+ZCOSSTRETCHLA*ZCOSLA*ZX) + ! (1-c.c)+(1+c.c)((sin lp.sin l + cos lp.cos l.cos(Lp-L)) + ! sin lr = ------------------------------------------------------- + ! D + ZSINLAS = (ZB + ZA*(ZSINSTRETCHLA*ZSINLA+ZCOSSTRETCHLA*ZCOSLA*ZX)) / ZD + ! D' = D * cos lr + ZD = ZD * (AMAX1(1e-6,SQRT(1.-ZSINLAS*ZSINLAS))) + ! 2.c.(cos lp.sin l - sin lp.cos l.cos(Lp-L)) + ! cos Lr = ------------------------------------------- + ! D' + ZCOSLOS = 2.*PCOEF*(ZCOSSTRETCHLA*ZSINLA-ZSINSTRETCHLA*ZCOSLA*ZX) / ZD + ! 2.c.cos l.cos(Lp-L) + ! sin Lr = ------------------- + ! D' + ZSINLOS = 2.*PCOEF*(ZCOSLA*ZY) / ZD + ! saturations (corrects calculation errors) + ZSINLAS = MAX(ZSINLAS,-1.) + ZSINLAS = MIN(ZSINLAS, 1.) + ZCOSLOS = MAX(ZCOSLOS,-1.) + ZCOSLOS = MIN(ZCOSLOS, 1.) + ! back from sine & cosine + PLAC(JLOOP1) = ASIN(ZSINLAS) + IF (ZSINLOS>0) THEN + PLOC(JLOOP1) = ACOS(ZCOSLOS) + ELSE + PLOC(JLOOP1) = -ACOS(ZCOSLOS) + ENDIF + ENDDO + PLOC(:) = PLOC(:) * 180. / XPI + PLAC(:) = PLAC(:) * 180. / XPI +IF (LHOOK) CALL DR_HOOK('ARPEGE_STRETCH_A',1,ZHOOK_HANDLE) +END SUBROUTINE ARPEGE_STRETCH_A +!------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_GAUSS_1COV diff --git a/src/SURFEX/hor_interpol_latlon_1cov.F90 b/src/SURFEX/hor_interpol_latlon_1cov.F90 new file mode 100644 index 000000000..d956e0517 --- /dev/null +++ b/src/SURFEX/hor_interpol_latlon_1cov.F90 @@ -0,0 +1,91 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_LATLON_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_LATLON_1COV* - Interpolation from a lat/lon regular grid on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_LATLON) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT, LINTERP +USE MODD_GRID_LATLONREGUL, ONLY : XILAT1, XILON1, XILAT2, XILON2,& + NINLAT, NINLON, NILENGTH,XILATARRAY +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODI_ADAPT_HORIBL_SURF +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKIN ! input mask +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKOUT ! output mask +INTEGER :: INO ! output number of points +INTEGER :: JL ! loop counter +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------------- +!* 1. Allocations +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_LATLON_1COV',0,ZHOOK_HANDLE) +INO = SIZE(XLAT_OUT) +! +ALLOCATE(IMASKIN (NILENGTH)) +! +ALLOCATE(IMASKOUT(INO)) +! +!* 2. Initializations +! +IMASKOUT = 1 +! +! +!* 3. Input mask +! +IMASKIN(:) = 1 +WHERE(PFIELDIN(:)==XUNDEF) IMASKIN = 0 +! +! +!* 4. Interpolation with horibl +! +CALL ADAPT_HORIBL_SURF(XILATARRAY,XILAT1,XILON1,XILAT2,XILON2,NINLAT,NINLON,NILENGTH, & + PFIELDIN(:),INO,XLON_OUT,XLAT_OUT,PFIELDOUT(:),.FALSE., & + KLUOUT,LINTERP,IMASKIN,IMASKOUT) +! +!* 6. Deallocations +! +IF (ALLOCATED(IMASKIN )) DEALLOCATE(IMASKIN ) +IF (ALLOCATED(IMASKOUT)) DEALLOCATE(IMASKOUT) +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_LATLON_1COV',1,ZHOOK_HANDLE) +!------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_LATLON_1COV diff --git a/src/SURFEX/hor_interpol_none_1cov.F90 b/src/SURFEX/hor_interpol_none_1cov.F90 new file mode 100644 index 000000000..27d3245b9 --- /dev/null +++ b/src/SURFEX/hor_interpol_none_1cov.F90 @@ -0,0 +1,78 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_NONE_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_NONE_1COV * - Only extrapolation on 1 vertical level only +!! +!! PURPOSE +!! ------- +! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_NONE) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT, CMASK +! +USE MODD_GRID_BUFFER, ONLY : NNI +USE MODD_SURF_PAR, ONLY : XUNDEF +! +USE MODI_PACK_SAME_RANK +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +USE MODI_GET_SURF_MASK_n +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, POINTER, DIMENSION(:) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! +INTEGER, DIMENSION(:), ALLOCATABLE :: IMASKOUT ! output mask +INTEGER :: INO ! output number of points +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------ +! +!* 1. Initialisation of the output mask +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_NONE_1COV',0,ZHOOK_HANDLE) +INO = SIZE(XLAT_OUT) +ALLOCATE(IMASKOUT(INO)) + CALL GET_SURF_MASK_n(CMASK,INO,IMASKOUT,NNI,KLUOUT) +! +!* 2. Mask the input field with the output mask +!!mask du tableau de taille FULL en fonction du type de surface + CALL PACK_SAME_RANK(IMASKOUT,PFIELDIN,PFIELDOUT) +! +!* 6. Deallocations +! +DEALLOCATE(IMASKOUT) +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_NONE_1COV',1,ZHOOK_HANDLE) + +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_NONE_1COV diff --git a/src/SURFEX/hor_interpol_rotlatlon_1cov.F90 b/src/SURFEX/hor_interpol_rotlatlon_1cov.F90 new file mode 100644 index 000000000..32ea20a62 --- /dev/null +++ b/src/SURFEX/hor_interpol_rotlatlon_1cov.F90 @@ -0,0 +1,305 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ######### +SUBROUTINE HOR_INTERPOL_ROTLATLON_1COV(KLUOUT,PFIELDIN,PFIELDOUT) +! ################################################################################# +! +!!**** *HOR_INTERPOL_ROTLATLON_1COV * - Interpolation from a rotated lat/lon grid on 1 vertical level only +!! +!! PURPOSE +!! ------- +!! +!!** METHOD +!! ------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M.Moge (adapted from HOR_INTERPOL_ROTLATLON) +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/2015 +!!------------------------------------------------------------------ +! +! +! +USE MODD_PREP, ONLY : XLAT_OUT, XLON_OUT +USE MODD_GRID_ROTLATLON +USE MODD_SURF_PAR, ONLY : XUNDEF +USE MODD_GRID_GRIB, ONLY : NNI + +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! +INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing +REAL, DIMENSION(:), INTENT(IN) :: PFIELDIN ! field to interpolate horizontally +REAL, DIMENSION(:), INTENT(OUT) :: PFIELDOUT ! interpolated field +! +!* 0.2 declarations of local variables +! + + INTEGER, ALLOCATABLE :: ii(:),jj(:) + + REAL, ALLOCATABLE :: XLAT_IND(:),XLON_IND(:), & + XRAT_OUT(:),XRON_OUT(:), & + w00(:),w01(:), & + w10(:),w11(:) + + LOGICAL, ALLOCATABLE :: LMASK(:) + +INTEGER :: I,J,K,L,IJ,IJ00,IJ01,IJ10,IJ11,INO,JL +REAL :: WX,WY,WSUM +REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!------------------------------------------------------------------------------------- +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_ROTLATLON_1COV',0,ZHOOK_HANDLE) +WRITE(KLUOUT,'(A)')' | Running rotated latlon interpolation' + +INO = SIZE(XLAT_OUT) + +! +!* 1. Allocations +! +ALLOCATE(XRAT_OUT(INO), & + XRON_OUT(INO), & + XLAT_IND(INO), & + XLON_IND(INO), & + II(INO), & + JJ(INO), & + W00(INO), & + W01(INO), & + W10(INO), & + W11(INO)) + +ALLOCATE(LMASK(NNI)) +! +!* Transformation of latitudes/longitudes into rotated coordinates + + WRITE(KLUOUT,*)'XLAT_OUT',XLAT_OUT(10:10) + WRITE(KLUOUT,*)'XLON_OUT',XLON_OUT(10:10) + + CALL REGROT(XLON_OUT,XLAT_OUT, & + XRON_OUT,XRAT_OUT, & + INO,1,INO,1, & + XRLOP,XRLAP,1) + + WRITE(KLUOUT,*)'XRAT_OUT',XRAT_OUT(10:10) + WRITE(KLUOUT,*)'XRON_OUT',XRON_OUT(10:10) + + DO IJ=1,INO + XLAT_IND(IJ) = ( XRAT_OUT(IJ) - XRILA1) / XRDY + 1. + XLON_IND(IJ) = ( XRON_OUT(IJ) - XRILO1) / XRDX + 1. + ENDDO + + PFIELDOUT(:) = XUNDEF + + + LMASK= .TRUE. + WHERE ( ABS(PFIELDIN(:)-XUNDEF) < 1.e-6 ) LMASK = .FALSE. + + DO IJ=1,INO + + II(IJ) = INT(XLON_IND(IJ)) + JJ(IJ) = INT(XLAT_IND(IJ)) + + WX = XLON_IND(IJ) - FLOAT(II(IJ)) + WY = XLAT_IND(IJ) - FLOAT(JJ(IJ)) + + W00(IJ) = (1.-WX)*(1.-WY) + W01(IJ) = (1.-WX)* WY + W10(IJ) = WX *(1.-WY) + W11(IJ) = WX * WY + + K = II(IJ) + L = JJ(IJ) + IJ00 = k + NRX*(l -1) + IJ01 = k + NRX*(l+1 -1) + IJ10 = k+1 + NRX*(l -1) + IJ11 = k+1 + NRX*(l+1 -1) + + IF (.NOT. LMASK(IJ00)) w00(IJ) = 0. + IF (.NOT. LMASK(IJ01)) w01(IJ) = 0. + IF (.NOT. LMASK(IJ10)) w10(IJ) = 0. + IF (.NOT. LMASK(IJ11)) w11(IJ) = 0. + + wsum = w00(IJ) + w01(IJ) + & + w10(IJ) + w11(IJ) + + IF ( ABS(wsum) < 1.e-6 ) CYCLE + + w00(IJ) = w00(IJ) / wsum + w01(IJ) = w01(IJ) / wsum + w10(IJ) = w10(IJ) / wsum + w11(IJ) = w11(IJ) / wsum + + ENDDO + + ! + ! Bi linear + ! + + WRITE(KLUOUT,*)'NRX,NRY',NRX,NRY + + DO IJ=1,INO + + K = II(IJ) + L = JJ(IJ) + IJ00 = k + NRX*(l -1) + IJ01 = k + NRX*(l+1 -1) + IJ10 = k+1 + NRX*(l -1) + IJ11 = k+1 + NRX*(l+1 -1) + + WRITE(KLUOUT,*)PFIELDIN(IJ00) + + PFIELDOUT(IJ) = w00(IJ)*PFIELDIN(IJ00) + & + w01(IJ)*PFIELDIN(IJ01) + & + w10(IJ)*PFIELDIN(IJ10) + & + w11(IJ)*PFIELDIN(IJ11) + + ENDDO + + +! +!* 5. Deallocations +! +DEALLOCATE(XRAT_OUT,XRON_OUT, & + XLAT_IND,XLON_IND, & + II,JJ, & + W00,W01,W10,W11, & + LMASK) +! +IF (LHOOK) CALL DR_HOOK('HOR_INTERPOL_ROTLATLON_1COV',1,ZHOOK_HANDLE) +CONTAINS +! + SUBROUTINE REGROT(PXREG,PYREG,PXROT,PYROT,KXDIM,KYDIM,KX,KY, & + PXCEN,PYCEN,KCALL) +! +! + IMPLICIT NONE +! +!----------------------------------------------------------------------- +! +!* CONVERSION BETWEEN REGULAR AND ROTATED SPHERICAL COORDINATES. +!* +!* PXREG LONGITUDES OF THE REGULAR COORDINATES +!* PYREG LATITUDES OF THE REGULAR COORDINATES +!* PXROT LONGITUDES OF THE ROTATED COORDINATES +!* PYROT LATITUDES OF THE ROTATED COORDINATES +!* ALL COORDINATES GIVEN IN DEGREES N (NEGATIVE FOR S) +!* AND DEGREES E (NEGATIVE VALUES FOR W) +!* KXDIM DIMENSION OF THE GRIDPOINT FIELDS IN THE X-DIRECTION +!* KYDIM DIMENSION OF THE GRIDPOINT FIELDS IN THE Y-DIRECTION +!* KX NUMBER OF GRIDPOINT IN THE X-DIRECTION +!* KY NUMBER OF GRIDPOINTS IN THE Y-DIRECTION +!* PXCEN REGULAR LONGITUDE OF THE SOUTH POLE OF THE ROTATED GRID +!* PYCEN REGULAR LATITUDE OF THE SOUTH POLE OF THE ROTATED GRID +!* +!* KCALL=-1: FIND REGULAR AS FUNCTIONS OF ROTATED COORDINATES. +!* KCALL= 1: FIND ROTATED AS FUNCTIONS OF REGULAR COORDINATES. +!* +!* J.E. HAUGEN HIRLAM JUNE -92 +! +!----------------------------------------------------------------------- +! + INTEGER KXDIM,KYDIM,KX,KY,KCALL + REAL PXREG(KXDIM,KYDIM),PYREG(KXDIM,KYDIM), & + PXROT(KXDIM,KYDIM),PYROT(KXDIM,KYDIM), & + PXCEN,PYCEN +! +!----------------------------------------------------------------------- +! + REAL PI,ZRAD,ZSYCEN,ZCYCEN,ZXMXC,ZSXMXC,ZCXMXC,ZSYREG,ZCYREG, & + ZSYROT,ZCYROT,ZCXROT,ZSXROT,ZRADI + INTEGER JY,JX + REAL(KIND=JPRB) :: ZHOOK_HANDLE +! +!----------------------------------------------------------------------- +! + IF (LHOOK) CALL DR_HOOK('REGROT',0,ZHOOK_HANDLE) + PI = 4.*ATAN(1.) + ZRAD = PI/180. + ZRADI = 1./ZRAD + ZSYCEN = SIN(ZRAD*(PYCEN+90.)) + ZCYCEN = COS(ZRAD*(PYCEN+90.)) +! + IF (KCALL.EQ.1) THEN +! + DO JY = 1,KY + DO JX = 1,KX +! + ZXMXC = ZRAD*(PXREG(JX,JY) - PXCEN) + ZSXMXC = SIN(ZXMXC) + ZCXMXC = COS(ZXMXC) + ZSYREG = SIN(ZRAD*PYREG(JX,JY)) + ZCYREG = COS(ZRAD*PYREG(JX,JY)) + ZSYROT = ZCYCEN*ZSYREG - ZSYCEN*ZCYREG*ZCXMXC + ZSYROT = MAX(ZSYROT,-1.0) + ZSYROT = MIN(ZSYROT,+1.0) +! + PYROT(JX,JY) = ASIN(ZSYROT)*ZRADI +! + ZCYROT = COS(PYROT(JX,JY)*ZRAD) + ZCXROT = (ZCYCEN*ZCYREG*ZCXMXC + & + ZSYCEN*ZSYREG)/ZCYROT + ZCXROT = MAX(ZCXROT,-1.0) + ZCXROT = MIN(ZCXROT,+1.0) + ZSXROT = ZCYREG*ZSXMXC/ZCYROT +! + PXROT(JX,JY) = ACOS(ZCXROT)*ZRADI +! + IF (ZSXROT.LT.0.0) PXROT(JX,JY) = -PXROT(JX,JY) +! + ENDDO + ENDDO +! + ELSEIF (KCALL.EQ.-1) THEN +! + DO JY = 1,KY + DO JX = 1,KX +! + ZSXROT = SIN(ZRAD*PXROT(JX,JY)) + ZCXROT = COS(ZRAD*PXROT(JX,JY)) + ZSYROT = SIN(ZRAD*PYROT(JX,JY)) + ZCYROT = COS(ZRAD*PYROT(JX,JY)) + ZSYREG = ZCYCEN*ZSYROT + ZSYCEN*ZCYROT*ZCXROT + ZSYREG = MAX(ZSYREG,-1.0) + ZSYREG = MIN(ZSYREG,+1.0) +! + PYREG(JX,JY) = ASIN(ZSYREG)*ZRADI +! + ZCYREG = COS(PYREG(JX,JY)*ZRAD) + ZCXMXC = (ZCYCEN*ZCYROT*ZCXROT - & + ZSYCEN*ZSYROT)/ZCYREG + ZCXMXC = MAX(ZCXMXC,-1.0) + ZCXMXC = MIN(ZCXMXC,+1.0) + ZSXMXC = ZCYROT*ZSXROT/ZCYREG + ZXMXC = ACOS(ZCXMXC)*ZRADI + IF (ZSXMXC.LT.0.0) ZXMXC = -ZXMXC +! + PXREG(JX,JY) = ZXMXC + PXCEN +! + ENDDO + ENDDO +! + ELSE + WRITE(6,'(1X,''INVALID KCALL IN REGROT'')') + CALL ABORT + ENDIF + IF (LHOOK) CALL DR_HOOK('REGROT',1,ZHOOK_HANDLE) +! + END SUBROUTINE REGROT +! +!------------------------------------------------------------------------------------- +END SUBROUTINE HOR_INTERPOL_ROTLATLON_1COV diff --git a/src/SURFEX/mode_extend_grid_parameter.F90 b/src/SURFEX/mode_extend_grid_parameter.F90 new file mode 100644 index 000000000..fb81bf991 --- /dev/null +++ b/src/SURFEX/mode_extend_grid_parameter.F90 @@ -0,0 +1,117 @@ +!################## +MODULE MODE_EXTEND_GRID_PARAMETER +!################## +! +CONTAINS +! +! Author +! M.Moge 01/03/2015 +! +! ############################################################# + SUBROUTINE EXTEND_GRID_PARAMETERX1(HPROGRAM,HGRID,HREC,KDIM,KSIZE,KIMAX_ll,KJMAX_ll,PFIELD,PFIELD_EXTEND) +! ############################################################# +! +!!**** * - routine to extend a real splitted array on SURFEX halo +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +#ifdef OL +!USE MODE_EXTEND_GRID_PARAMETER_OL +#endif +#ifdef MNH +USE MODI_EXTEND_GRID_PARAMETERX1_MNH +#endif +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program + CHARACTER(LEN=10), INTENT(IN) :: HGRID ! grid type + CHARACTER(LEN=6), INTENT(IN) :: HREC ! name of the parameter +INTEGER, INTENT(IN) :: KDIM ! size of PFIELD +INTEGER, INTENT(IN) :: KSIZE ! size of PFIELD_EXTEND +INTEGER, INTENT(IN) :: KIMAX_ll !(global) dimension of the domain - X direction +INTEGER, INTENT(IN) :: KJMAX_ll !(global) dimension of the domain - Y direction +REAL, DIMENSION(KDIM ), INTENT(IN) :: PFIELD ! real field for complete grid +REAL, DIMENSION(KSIZE), INTENT(OUT):: PFIELD_EXTEND! real field for splitted grid +! +!* 0.2 Declarations of local variables +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE +!------------------------------------------------------------------------------- +IF (LHOOK) CALL DR_HOOK('MODE_EXTEND_GRID_PARAMETER:EXTEND_GRID_PARAMETERX1',0,ZHOOK_HANDLE) +! +IF (HPROGRAM=='MESONH') THEN +#ifdef MNH + CALL EXTEND_GRID_PARAMETERX1_MNH(HGRID,HREC,KDIM,KSIZE,KIMAX_ll,KJMAX_ll,PFIELD,PFIELD_EXTEND) +#endif +ENDIF +! +! +IF (HPROGRAM=='OFFLIN') THEN +#ifdef OL +! CALL EXTEND_GRID_PARAMETERX1_OL(HPROGRAM,HGRID,HREC,KDIM,KSIZE,PFIELD,PFIELD_EXTEND) +!TODO : write subroutine EXTEND_GRID_PARAMETERX1_OL +#endif +ENDIF +! +IF (LHOOK) CALL DR_HOOK('MODE_EXTEND_GRID_PARAMETER:EXTEND_GRID_PARAMETERX1',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +END SUBROUTINE EXTEND_GRID_PARAMETERX1 +! +! +! ############################################################# + SUBROUTINE EXTEND_GRID_PARAMETERN0(HPROGRAM,HGRID,HREC,KFIELD,KFIELD_EXTEND) +! ############################################################# +! +!!**** * - routine to "extend" a integer related to splitted grid on SURFEX halo +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +#ifdef OL +!USE MODE_EXTEND_GRID_PARAMETER_OL +#endif +#ifdef MNH +USE MODI_EXTEND_GRID_PARAMETERN0_MNH +#endif +! +IMPLICIT NONE +! +!* 0.1 Declarations of arguments +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program + CHARACTER(LEN=10), INTENT(IN) :: HGRID ! grid type + CHARACTER(LEN=6), INTENT(IN) :: HREC ! name of the parameter +INTEGER, INTENT(IN) :: KFIELD ! integer scalar for complete grid +INTEGER, INTENT(OUT):: KFIELD_EXTEND ! integer scalar for splitted grid +!* 0.2 Declarations of local variables +! +REAL(KIND=JPRB) :: ZHOOK_HANDLE +!------------------------------------------------------------------------------- +IF (LHOOK) CALL DR_HOOK('MODE_EXTEND_GRID_PARAMETER:EXTEND_GRID_PARAMETERN0',0,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +! +IF (HPROGRAM=='MESONH') THEN +#ifdef MNH + CALL EXTEND_GRID_PARAMETERN0_MNH(HGRID,HREC,KFIELD,KFIELD_EXTEND) +#endif +ENDIF +! +IF (HPROGRAM=='OFFLIN') THEN +#ifdef OL +! CALL EXTEND_GRID_PARAMETERN0_OL(HPROGRAM,HGRID,HREC,KFIELD,KFIELD_EXTEND) +!TODO : write subroutine EXTEND_GRID_PARAMETERN0_OL +#endif +ENDIF +! +IF (LHOOK) CALL DR_HOOK('MODE_EXTEND_GRID_PARAMETER:EXTEND_GRID_PARAMETERN0',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +END SUBROUTINE EXTEND_GRID_PARAMETERN0 +! +END MODULE MODE_EXTEND_GRID_PARAMETER diff --git a/src/SURFEX/read_covers_and_av_pgd_on_layers.F90 b/src/SURFEX/read_covers_and_av_pgd_on_layers.F90 new file mode 100644 index 000000000..220e086ee --- /dev/null +++ b/src/SURFEX/read_covers_and_av_pgd_on_layers.F90 @@ -0,0 +1,210 @@ +!SURFEX_LIC Copyright 1994-2014 Meteo-France +!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence +!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!SURFEX_LIC for details. version 1. +! ################################################################ + SUBROUTINE READ_COVERS_AND_AV_PGD_1D_ON_LAYERS(HFILEPGDTYPE,HRECFM,KLU,KDATA_LAYER,PFIELD2D,PDATA,HSFTYPE,HATYPE,PDZ,KDECADE) +! ################################################################ +!! +!! PURPOSE +!! ------- +!! +!! METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! M.Moge CNRS - LA +!! inspired from AV_PGD_1D +!! +!! MODIFICATION +!! ------------ +! +! +!! Original 06/05/2015 +!! +!---------------------------------------------------------------------------- +! +!* 0. DECLARATION +! ----------- +! +USE MODD_SURF_PAR, ONLY : XUNDEF +USE MODD_DATA_COVER, ONLY : XDATA_BLD_HEIGHT +USE MODD_DATA_COVER_n, ONLY : XDATA_NATURE, XDATA_TOWN, XDATA_BLD, XDATA_GARDEN, & + XDATA_SEA, XDATA_WATER, XDATA_VEGTYPE +USE MODD_DATA_COVER_PAR, ONLY : NVT_TREE, NVT_CONI, NVT_EVER, XCDREF, JPCOVER +! +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +USE MODI_ABOR1_SFX +! +USE MODI_READ_SURF +#ifdef MNH +USE MODI_READ_SURFX2COV_1COV_MNH +#endif +! +IMPLICIT NONE +! +!* 0.1 Declaration of arguments +! ------------------------ +! +CHARACTER(LEN=6), INTENT(IN) :: HFILEPGDTYPE ! type of input file +CHARACTER(LEN=12), INTENT(IN) :: HRECFM ! Name of the article to be read +INTEGER, INTENT(IN) :: KLU ! number of points +INTEGER, INTENT(IN) :: KDATA_LAYER ! number of layers +REAL, DIMENSION(KLU,KDATA_LAYER), INTENT(OUT) :: PFIELD2D ! secondary field to construct +REAL, DIMENSION(JPCOVER,KDATA_LAYER), INTENT(IN) :: PDATA ! secondary field value for each class +CHARACTER(LEN=3), INTENT(IN) :: HSFTYPE ! Type of surface where the field is defined +CHARACTER(LEN=3), INTENT(IN) :: HATYPE ! Type of averaging +REAL, DIMENSION(KLU), INTENT(IN), OPTIONAL :: PDZ ! first model half level +INTEGER, INTENT(IN), OPTIONAL :: KDECADE ! current month +! +!* 0.2 Declaration of local variables +! ------------------------------ +! +! +INTEGER :: ICOVER ! number of cover classes +INTEGER :: JCOVER ! loop on cover classes +INTEGER :: JLAYER ! loop on layers +! +REAL, DIMENSION(KLU) :: ZWORK, ZDZ +REAL :: ZWEIGHT +REAL, DIMENSION(KLU) :: ZCOVER_WEIGHT +REAL :: ZDATA_COVER +REAL, DIMENSION(KLU) :: ZSUM_COVER_WEIGHT +REAL, DIMENSION(KLU) :: ZWEIGHT_MAX +REAL(KIND=JPRB) :: ZHOOK_HANDLE +LOGICAL, DIMENSION(JPCOVER) :: GCOVER ! flag to read the covers +REAL, DIMENSION(KLU) :: ZCOVER ! cover fractions +CHARACTER(LEN=100) :: YCOMMENT +INTEGER :: IRESP ! reading return code +CHARACTER(LEN=16) :: YRECFM ! Name of the article to be read +!------------------------------------------------------------------------------- +! +IF (LHOOK) CALL DR_HOOK('READ_COVERS_AND_AV_PGD_1D_ON_LAYERS',0,ZHOOK_HANDLE) +! +!* 0.3 Initializations +! +IF (PRESENT(PDZ)) THEN + ZDZ(:)=PDZ(:) +ELSE + ZDZ(:)=XCDREF +END IF +! +PFIELD2D(:,:)=XUNDEF +! +! +!* depths are deduced from the cover types +!* reading of the cover to obtain the thickness of layers +CALL READ_SURF(HFILEPGDTYPE,HRECFM,GCOVER(:),IRESP,HDIR='-') +YRECFM='COVER' +#ifdef MNH +! +! Loop on layers +DO JLAYER=1,KDATA_LAYER + ZWORK(:)=0. + ZWEIGHT_MAX(:)=0. + ZSUM_COVER_WEIGHT(:)=0. + ! loop on covers + DO JCOVER=1,JPCOVER + ! + !* 1. depths are deduced from the cover types + ! reading of the cover to obtain the thickness of layers + ! + IF ( GCOVER( JCOVER ) ) THEN + CALL READ_SURFX2COV_1COV_MNH(YRECFM,KLU,JCOVER,ZCOVER(:),IRESP,YCOMMENT,'A') + ELSE + ZCOVER(:) = 0. + ENDIF + ! + !* 2. averaging + ! + ! 2.1. Selection of the weighting function + SELECT CASE (HSFTYPE) + CASE('ALL') + ZWEIGHT=1. + CASE('NAT') + ZWEIGHT=XDATA_NATURE(JCOVER) + CASE('GRD') + ZWEIGHT=XDATA_TOWN (JCOVER) * XDATA_GARDEN(JCOVER) + CASE('TWN') + ZWEIGHT=XDATA_TOWN (JCOVER) + CASE('WAT') + ZWEIGHT=XDATA_WATER (JCOVER) + CASE('SEA') + ZWEIGHT=XDATA_SEA (JCOVER) + CASE('BLD') + ZWEIGHT=XDATA_TOWN (JCOVER) * XDATA_BLD(JCOVER) + CASE('BLV') !* building Volume + ZWEIGHT=XDATA_TOWN (JCOVER) * XDATA_BLD(JCOVER) * XDATA_BLD_HEIGHT(JCOVER) + CASE('STR') + ZWEIGHT=XDATA_TOWN (JCOVER) * ( 1. - XDATA_BLD(JCOVER) ) + CASE('TRE') + ZWEIGHT=XDATA_NATURE(JCOVER) * ( XDATA_VEGTYPE(JCOVER,NVT_TREE) + XDATA_VEGTYPE(JCOVER,NVT_EVER) + XDATA_VEGTYPE(JCOVER,NVT_CONI) ) + CASE('GRT') + ZWEIGHT=XDATA_TOWN(JCOVER) * XDATA_GARDEN(JCOVER) * ( XDATA_VEGTYPE(JCOVER,NVT_TREE) + XDATA_VEGTYPE(JCOVER,NVT_EVER) + XDATA_VEGTYPE(JCOVER,NVT_CONI) ) + CASE DEFAULT + CALL ABOR1_SFX('AV_PGD_1D: WEIGHTING FUNCTION NOT ALLOWED '//HSFTYPE) + END SELECT + ! 2.2. Averaging + ZCOVER_WEIGHT(:) = ZCOVER(:) * ZWEIGHT + ZSUM_COVER_WEIGHT(:) = ZSUM_COVER_WEIGHT(:) + ZCOVER_WEIGHT(:) + ZDATA_COVER = PDATA(JCOVER,JLAYER) + SELECT CASE (HATYPE) + CASE ('ARI') + ZWORK(:) = ZWORK(:) + ZDATA_COVER * ZCOVER_WEIGHT(:) + CASE('INV' ) + ZWORK (:)= ZWORK(:) + 1./ZDATA_COVER * ZCOVER_WEIGHT(:) + CASE('CDN') + ZWORK (:)= ZWORK(:) + 1./(LOG(ZDZ(:)/ZDATA_COVER))**2 * ZCOVER_WEIGHT(:) + CASE('MAJ' ) + WHERE(ZCOVER_WEIGHT(:)>ZWEIGHT_MAX(:)) + ZWEIGHT_MAX(:) = ZCOVER_WEIGHT(:) + ZWORK (:) = ZDATA_COVER + END WHERE + CASE DEFAULT + CALL ABOR1_SFX('AV_PGD_1D: (1) AVERAGING TYPE NOT ALLOWED : "'//HATYPE//'"') + END SELECT + ! + END DO ! DO JCOVER=1,JPCOVER + ! + ! 2.3. End of Averaging + SELECT CASE (HATYPE) + CASE ('ARI') + WHERE ( ZSUM_COVER_WEIGHT(:) >0. ) + PFIELD2D(:,JLAYER) = ZWORK(:) / ZSUM_COVER_WEIGHT(:) + END WHERE + CASE('INV' ) + WHERE ( ZSUM_COVER_WEIGHT(:) >0. ) + PFIELD2D(:,JLAYER) = ZSUM_COVER_WEIGHT(:) / ZWORK(:) + END WHERE + CASE('CDN') + WHERE ( ZSUM_COVER_WEIGHT(:) >0. ) + PFIELD2D(:,JLAYER) = ZDZ(:) * EXP( - SQRT(ZSUM_COVER_WEIGHT(:)/ZWORK(:)) ) + END WHERE + CASE('MAJ' ) + WHERE ( ZSUM_COVER_WEIGHT(:) >0. ) + PFIELD2D(:,JLAYER) = ZWORK(:) + END WHERE + CASE DEFAULT + CALL ABOR1_SFX('AV_PGD_1D: (2) AVERAGING TYPE NOT ALLOWED') + END SELECT +! +END DO !DO JLAYER=1,KDATA_LAYER +#endif +! +IF (LHOOK) CALL DR_HOOK('READ_COVERS_AND_AV_PGD_1D_ON_LAYERS',1,ZHOOK_HANDLE) +!------------------------------------------------------------------------------- +END SUBROUTINE READ_COVERS_AND_AV_PGD_1D_ON_LAYERS diff --git a/src/SURFEX/write_lcover.F90 b/src/SURFEX/write_lcover.F90 new file mode 100644 index 000000000..853c8a89a --- /dev/null +++ b/src/SURFEX/write_lcover.F90 @@ -0,0 +1,86 @@ +! ######### + SUBROUTINE WRITE_LCOVER(HPROGRAM,OCOVER) +! ################################ +! +!!**** *READ_LCOVER* - routine to write a file for +!! physiographic data file of model _n +!! +!! PURPOSE +!! ------- +!! The purpose of this routine is to write the list of covers to a file in parallel using MPI +!! +!! +!!** METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! M. Moge *LA - CNRS* +!! +!! MODIFICATIONS +!! ------------- +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_DATA_COVER_PAR, ONLY : JPCOVER +!USE MODD_WATFLUX_n, ONLY : LCOVER +! +USE MODI_WRITE_SURF +! +USE YOMHOOK ,ONLY : LHOOK, DR_HOOK +USE PARKIND1 ,ONLY : JPRB +! +IMPLICIT NONE +! +#ifndef NOMPI +INCLUDE "mpif.h" +#endif +! +!* 0.1 Declarations of arguments +! ------------------------- +! + CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program +LOGICAL, DIMENSION(JPCOVER) :: OCOVER ! list of covers +! +!* 0.2 Declarations of local variables +! ------------------------------- +! +INTEGER :: IRESP ! Error code after reading +CHARACTER(LEN=12) :: YRECFM ! Name of the article to be read +CHARACTER(LEN=100):: YCOMMENT ! Comment string +LOGICAL, DIMENSION(JPCOVER) :: GCOVER ! tmp list of covers +REAL(KIND=JPRB) :: ZHOOK_HANDLE +INTEGER :: IINFO +!------------------------------------------------------------------------------- +! +! +!* ascendant compatibility +IF (LHOOK) CALL DR_HOOK('WRITE_LCOVER',0,ZHOOK_HANDLE) +#ifndef NOMPI +CALL MPI_ALLREDUCE(OCOVER, GCOVER, SIZE(OCOVER),MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, IINFO) +#endif +OCOVER(:)=GCOVER(:) +YRECFM='COVER_LIST' +YCOMMENT='(LOGICAL LIST)' +CALL WRITE_SURF(HPROGRAM,YRECFM,OCOVER(:),IRESP,HCOMMENT=YCOMMENT,HDIR='-') +! +IF (LHOOK) CALL DR_HOOK('WRITE_LCOVER',1,ZHOOK_HANDLE) +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE WRITE_LCOVER -- GitLab