!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: /srv/cvsroot/MNH-VX-Y-Z/src/MNH/richardson.f90,v $ $Revision: 1.1.10.1.18.2 $ ! MASDEV4_7 solver 2006/05/18 13:07:25 !----------------------------------------------------------------- ! ######spl MODULE MODI_RICHARDSONZ2 ! ###################### ! INTERFACE ! SUBROUTINE RICHARDSONZ2(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ, & PRHODJ, PRHODREF,PTHETAV,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, & KIFAXX,KIFAXY,KITR,KTCOUNT,PY,PPHI,PBFB,& PBF_SXP2_YP1_Z) ! IMPLICIT NONE ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type ! ! Metric coefficients: REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! density of reference REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at ! time t ! REAL, INTENT(IN) :: PDXHATM ! mean grid increment in the x ! direction REAL, INTENT(IN) :: PDYHATM ! mean grid increment in the y ! direction ! REAL, DIMENSION (:), INTENT(IN) :: PRHOM ! mean of XRHODJ on the plane x y ! localized at a mass level ! REAL, DIMENSION(:), INTENT(IN) :: PAF,PCF ! vectors giving the nonvanishing REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF ! elements of the tri-diag. ! matrix in the pressure eq. ! ! arrays of sin or cos values ! for the FFT : REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX ! - along x REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY ! - along y ! ! decomposition in prime ! numbers for the FFT: INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX ! - along x INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY ! - along y ! INTEGER, INTENT(IN) :: KITR ! number of iterations for the ! pressure solver INTEGER, INTENT(IN) :: KTCOUNT ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation !JUAN Z_SPLITING REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB ! elements of the tri-diag. b-slide ! matrix in the pressure eq. REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide ! matrix in the pressure eq. !JUAN Z_SPLITING ! END SUBROUTINE RICHARDSONZ2 ! END INTERFACE ! END MODULE MODI_RICHARDSONZ2 ! ######spl SUBROUTINE RICHARDSONZ2(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ, & PRHODJ,PRHODREF,PTHETAV,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, & KIFAXX,KIFAXY,KITR,KTCOUNT,PY,PPHI,PBFB,& PBF_SXP2_YP1_Z) ! ######################################################################## ! !!**** *RICHARDSONZ2 * - solve the elliptic pressure equation by the Richardson's !! method !! !! PURPOSE !! ------- ! The purpose of this routine is to solve the elliptic pressure equation ! by the Richardson's method preconditioned by the flat Laplacian. ! !!** METHOD !! ------ !! !! The equation to be solved reads: !! !! Q (PHI) = Y !! !! where Q is the quasi-Laplacian ( subroutine QLAP ) and PHI the pressure !! function. !! We precondition the problem by the operator F : !! -1 -1 !! F * Q (PHI) = F (Y) !! F represents the flat quasi-laplacian ie. without orography. Its !! inversion is realized in the routine FLAT_INV. !! This equation is solved with the Richardson's method: !! !! (n+1) (n) ( -1 -1 (n) ) !! PHI = PHI + RELAX ( F ( Y ) - F * Q ( PHI ) ) !! ( ) !! !! EXTERNAL !! -------- !! Subroutine GDIV: compute J times the divergence of 1/J times a vector !! Function QLAP: compute the complete quasi-Laplacian Q !! Subroutine FLAT_INV : invert the flat quasi-laplacien F !! !! IMPLICIT ARGUMENTS !! ------------------ !! Module MODD_CONF: model configuration !! CCONF: option for the first time step !! of the segment ( start or restart) !! !! REFERENCE !! --------- !! Book2 of documentation (routine RICHARDSONZ2) !! !! AUTHOR !! ------ !! P. HÅreil and J. Stein * Meteo France * !! !! MODIFICATIONS !! ------------- !! Original 25/07/94 !! 14/01/97 Durran anelastic equation (Stein,Lafore) !! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS ! ------------ USE MODE_ll USE MODE_IO_ll USE MODD_CONF USE MODI_GDIV USE MODI_QLAP USE MODI_FLAT_INVZ USE MODD_VAR_ll, ONLY: IP USE MODD_IBM_PARAM_n USE MODI_SECOND_MNH ! IMPLICIT NONE ! !* 0.1 declarations of arguments ! ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type ! ! Metric coefficients: REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! density of reference REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at ! time t REAL, INTENT(IN) :: PDXHATM ! mean grid increment in the x ! direction REAL, INTENT(IN) :: PDYHATM ! mean grid increment in the y ! direction ! REAL, DIMENSION (:), INTENT(IN) :: PRHOM ! mean of XRHODJ on the plane x y ! localized at a mass level ! REAL, DIMENSION(:), INTENT(IN) :: PAF,PCF ! vectors giving the nonvanishing REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF ! elements of the tri-diag. ! matrix in the pressure eq. ! ! arrays of sin or cos values ! for the FFT : REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX ! - along x REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY ! - along y ! ! decomposition in prime ! numbers for the FFT: INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX ! - along x INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY ! - along y ! INTEGER, INTENT(IN) :: KITR ! number of iterations for the ! pressure solver INTEGER, INTENT(IN) :: KTCOUNT ! REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation ! !* 0.2 declarations of local variables ! INTEGER :: JM ! loop index ! REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZF_1_Y ! RHS of the preconditioned problem ! REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZCORREC ! iterative correction to the solution ! REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZWORK,ZRESIDUE ! quasi-laplacien Q of the iterative solution !JUAN Z_SPLITING REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB ! elements of the tri-diag. b-slide ! matrix in the pressure eq. REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide ! matrix in the pressure eq. !JUAN Z_SPLITING REAL :: ZMAX_RESIDUE,ZMAX_RESIDUE1,ZMAX_RESIDUE2,ZMAX_RESIVAR,ZMAX_RESIVAR1,ZMAX_RESIVAR2,ZMAX_Y INTEGER :: IINFO_ll INTEGER :: ILUIBMPRES,IRESPIBMPRES REAL :: ZTIME,ZTIME1,ZTIME2,ZTIME3,ZTIME4 ! ! !------------------------------------------------------------------------------- ! ILUIBMPRES = 103!FORT IF ((KTCOUNT==0.and.IP==1).OR.(IP==1.AND.KTCOUNT==1.AND.CCONF=='RESTA')) THEN OPEN(UNIT=ILUIBMPRES , FILE='ibm_rich.dat', IOSTAT=IRESPIBMPRES , FORM='FORMATTED' , & STATUS='NEW', POSITION='APPEND', ACCESS='SEQUENTIAL', ACTION='WRITE') ELSEIF (mod(KTCOUNT,20)==0.and.IP==1) THEN OPEN(UNIT=ILUIBMPRES , FILE='ibm_rich.dat', IOSTAT=IRESPIBMPRES , FORM='FORMATTED' , & STATUS='OLD', POSITION='APPEND', ACCESS='SEQUENTIAL', ACTION='WRITE') ENDIF CALL SECOND_MNH(ZTIME1) ZMAX_Y = MAX_ll(ABS(PY),IINFO_ll) ZMAX_Y = ZMAX_Y/XIBM_DIV ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY ZMAX_RESIDUE2 = MAX_ll(ABS(ZRESIDUE),IINFO_ll) ZMAX_RESIDUE = 1. ZMAX_RESIVAR = ZMAX_RESIDUE2/XIBM_DIV JM = 1 CALL SECOND_MNH(ZTIME2) IF (IP==1.and.mod(KTCOUNT,20)==0) WRITE(UNIT=ILUIBMPRES,FMT='(5(1X,E16.8))') KTCOUNT*1.,JM*1.,ZMAX_Y,ZMAX_RESIVAR,ZMAX_RESIDUE IF (ZMAX_Y.lt.XIBM_RES .or. KITR==0) RETURN ! CALL SECOND_MNH(ZTIME3) CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZCORREC,PBFB,& PBF_SXP2_YP1_Z) ! PPHI = PPHI - ZCORREC ! !------------------------------------------------------------------------------- ! !* 2. ITERATIVE LOOP ! -------------- ! CALL SECOND_MNH(ZTIME4) ZTIME=ZTIME4+ZTIME2-ZTIME3-ZTIME1 DO WHILE ((ZMAX_RESIVAR.gt.XIBM_RES) .and. (ZMAX_RESIDUE.gt.XIBM_RES).and. (JM