| 1 |
C $Id: pfind.F,v 1.1 1998/05/25 20:21:05 cnh Exp $ |
| 2 |
#include "CPP_OPTIONS.h" |
| 3 |
#include "CPP_MACROS.h" |
| 4 |
# define _D(a) |
| 5 |
C===================================================================== |
| 6 |
C||| Procedure name: PFIND ||| |
| 7 |
C||| Function: Diagnose new pressure field such that ||| |
| 8 |
C||| del^2(PS) = RHS2D - del^2(PH) ||| |
| 9 |
C||| Comments: ||| |
| 10 |
C===================================================================== |
| 11 |
CStartOfInterface |
| 12 |
SUBROUTINE PFIND( RHS2D, rPrmBz, divH, |
| 13 |
O PH , PS ) |
| 14 |
IMPLICIT NONE |
| 15 |
C /--------------------------------------------------------------\ |
| 16 |
C | ==== Global data ========== | |
| 17 |
C \--------------------------------------------------------------/ |
| 18 |
#include "SIZE.h" |
| 19 |
#include "AJAINF.h" |
| 20 |
#include "OPERATORS.h" |
| 21 |
#include "GRID.h" |
| 22 |
#include "PARAMS.h" |
| 23 |
#include "OLDG.h" |
| 24 |
#include "MASKS.h" |
| 25 |
#include "CG2DA.h" |
| 26 |
C /--------------------------------------------------------------\ |
| 27 |
C | ==== Routine arguments ==== | |
| 28 |
C | RHS2D - div(GuHat) - 1/dt*div(UHat) | |
| 29 |
C | rprmBz - Density anomaly wrt inert reference (kg/m**3) | |
| 30 |
C | divH - Vertically integrated divergent flow (m**3/s) | |
| 31 |
C | PH - Hydrostatic "pressure" (m). | |
| 32 |
C | PS - Surface "pressure" (m). | |
| 33 |
C \--------------------------------------------------------------/ |
| 34 |
REAL rPrmBz(Nx,Ny,Nz) |
| 35 |
REAL RHS2D (Nx,Ny ) |
| 36 |
REAL divH (Nx,Ny ) |
| 37 |
REAL PH (Nx,Ny,Nz) |
| 38 |
REAL PS (Nx,Ny ) |
| 39 |
CEndOfInterface |
| 40 |
C /--------------------------------------------------------------\ |
| 41 |
C | ==== Local variables ====== | |
| 42 |
C | PHxhat - Vertically integrated d/dx PH. | |
| 43 |
C | PHyhat - Vertically integrated d/dy PH. | |
| 44 |
C | divHMx - Maximum vert. intergrated divergence (s**-1). | |
| 45 |
C | MAXIT - Limit on number of CG iterations. | |
| 46 |
C | rMax - Max. absolute value of elliptic eqn. rhs. | |
| 47 |
C | toler - tolerance for CG convergence. | |
| 48 |
C | freqCheck - Interval at which CG convergence is tested. | |
| 49 |
C | rStart2 - Starting CG residual. | |
| 50 |
C | rFinal2 - Converged CG residual. | |
| 51 |
C | numIt2 - No. of CG iterations to converge. | |
| 52 |
C \--------------------------------------------------------------/ |
| 53 |
REAL PHxhat (Nx+1,Ny+1) |
| 54 |
REAL PHyhat (Nx+1,Ny+1) |
| 55 |
REAL divHMx |
| 56 |
INTEGER MAXIT |
| 57 |
REAL rMax |
| 58 |
REAL toler |
| 59 |
INTEGER freqCheck |
| 60 |
INTEGER numIt2 |
| 61 |
REAL rStart2 |
| 62 |
REAL rFinal2 |
| 63 |
C /--------------------------------------------------------------\ |
| 64 |
C | Loop counters | |
| 65 |
C \--------------------------------------------------------------/ |
| 66 |
INTEGER I, J, K |
| 67 |
|
| 68 |
C _D(( ' S/R PFIND: MAXVAL(RHS2d) = ',MAXVAL(RHS2d))) |
| 69 |
c WRITE(0,*) ' S/R PFIND: MAXVAL(RHS2d) = ',MAXVAL(RHS2d) |
| 70 |
|
| 71 |
DO J=1,Ny |
| 72 |
DO I=1,Nx |
| 73 |
C /-------------------------------------------------------------\ |
| 74 |
C | RHS2d <- RHS2d - (Dx{PHxHat}+Dy{PHyHat})/RONIL | |
| 75 |
C | + fsFac*ZA*PS/DT**2 | |
| 76 |
C \-------------------------------------------------------------/ |
| 77 |
RHS2d(I,J) = RHS2d(I,J) |
| 78 |
& -ZA(_I3(1,I,J))*freeSurfFac*PS(I,J)/DELT/DELT |
| 79 |
RHS2d(I,J) = RHS2d(I,J)*PMASK(I,J,1)/aNorm2d |
| 80 |
ENDDO |
| 81 |
ENDDO |
| 82 |
divHMx = MAXVAL(divH*rPvolHat) |
| 83 |
IF ( divHMx .GT. divHMax ) THEN |
| 84 |
IF ( toler2d .GT. 1.E-15) toler2d = toler2d * 0.9 |
| 85 |
ELSEIF ( divHMx .LT. divHMin ) THEN |
| 86 |
IF ( toler2d .LT. 1. ) toler2d = toler2d * 1.1 |
| 87 |
ENDIF |
| 88 |
c WRITE(0,*) ' S/R PFIND MAX(DIVxy) = ',divHMx |
| 89 |
c WRITE(0,*) ' S/R PFIND normalisation factor = ',1.D0/aNorm2d |
| 90 |
C RHS2d = RHS2d |
| 91 |
rMax = MAXVAL(ABS(RHS2d)) |
| 92 |
_D(( ' S/R PFIND: aNorm2d = ',aNorm2d)) |
| 93 |
_D(( ' S/R PFIND: rMax = ',rMax)) |
| 94 |
PSNM1 = PS |
| 95 |
IF ( rMax .EQ. 0. ) THEN |
| 96 |
PS = 0. |
| 97 |
ELSE |
| 98 |
RHS2d = RHS2d/rMax |
| 99 |
PS = PS/rMax |
| 100 |
MAXIT = MAX2DIT |
| 101 |
TOLER = TOLER2D |
| 102 |
FREQCHECK = MAX(freqCheckToler2d,1) |
| 103 |
C /--------------------------------------------------------------\ |
| 104 |
C | del2{PS} = RHS2d | |
| 105 |
C \--------------------------------------------------------------/ |
| 106 |
C OPEN (11,FILE='rhs',FORM='unformatted') |
| 107 |
C WRITE(11) RHS2d |
| 108 |
C CLOSE(11) |
| 109 |
|
| 110 |
CALL CG2D_OLD( |
| 111 |
C CALL CG2D ( |
| 112 |
I RHS2d, maxIt, freqCheck, toler, |
| 113 |
U PS, |
| 114 |
O numIt2, rStart2, rFinal2 |
| 115 |
& ) |
| 116 |
PS = PS*rMax |
| 117 |
IF ( numIt2 .GE. maxIt ) THEN |
| 118 |
WRITE(0,*) 'WARNING: 2d solver used maximum iterations.' |
| 119 |
WRITE(0,*) ' Requested tolerance = ',toler |
| 120 |
WRITE(0,*) ' Tolerance achieved = ',sqrt(rFinal2) |
| 121 |
WRITE(0,*) ' 2d iterations = ',numIt2 |
| 122 |
WRITE(0,*) ' Model timestep = ',nIter |
| 123 |
WRITE(0,*) ' Model time = ',currentTime |
| 124 |
ENDIF |
| 125 |
ENDIF |
| 126 |
_D(( ' S/R PFIND: MAXVAL(PS) = ',MAXVAL(PS) )) |
| 127 |
_D(( ' S/R PFIND: MINVAL(PS) = ',MINVAL(PS) )) |
| 128 |
CcnhDebugStarts |
| 129 |
c WRITE(0,*) ' S/R PFIND: nIts2, resid2 ', numIt2, sqrt(rFinal2) |
| 130 |
CcnhDebugEnds |
| 131 |
C |
| 132 |
RETURN |
| 133 |
END |
Parent Directory
| 
Revision Log
| 
Revision Graph