| 1 |
C $Id: update_t.F,v 1.10 1997/06/10 17:46:06 cnh Exp $ |
| 2 |
#include "CPP_OPTIONS.h" |
| 3 |
#include "CPP_MACROS.h" |
| 4 |
|
| 5 |
#define _DB0 0. |
| 6 |
|
| 7 |
C=================================================================== |
| 8 |
C Procedure name: UPDATE_T | |
| 9 |
C Function: Step forward the temperature field. | |
| 10 |
C Comments: | |
| 11 |
C=================================================================== |
| 12 |
CStartofinterface |
| 13 |
SUBROUTINE UPDATE_T( |
| 14 |
& U, V, W, K, T, GT ) |
| 15 |
IMPLICIT NONE |
| 16 |
C============== Global data ========================================== |
| 17 |
#include "SIZE.h" |
| 18 |
#include "AJAINF.h" |
| 19 |
#include "GRID.h" |
| 20 |
#include "PARAMS.h" |
| 21 |
#include "OPERATORS.h" |
| 22 |
#include "OLDG.h" |
| 23 |
#include "MASKS.h" |
| 24 |
#include "FORCING.h" |
| 25 |
C============= Routine arguments ===================================== |
| 26 |
C /--------------------------------------------------------------\ |
| 27 |
C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). | |
| 28 |
C | K - Working level. | |
| 29 |
C | T - Potential temperature (oC). | |
| 30 |
C \--------------------------------------------------------------/ |
| 31 |
REAL U (Nx,Ny,Nz) |
| 32 |
REAL V (Nx,Ny,Nz) |
| 33 |
REAL W (Nx,Ny,Nz) |
| 34 |
INTEGER K |
| 35 |
REAL T (Nx,Ny,Nz) |
| 36 |
REAL gt (Nx,Ny,Nz) |
| 37 |
CEndofinterface |
| 38 |
C============= Local variables ====================================== |
| 39 |
C /--------------------------------------------------------------\ |
| 40 |
C |Slice Notation | |
| 41 |
C |phiK - XY slice of variable "phi" at level K. | |
| 42 |
C |phiKP1 - XY slice of variable "phi" at level K+1. | |
| 43 |
C |phiKM1 - XY slice of variable "phi" at level K-1. | |
| 44 |
C |phiNM1 - Variable "phi" at time level N-1. | |
| 45 |
C |phiU - Varaible "phi" interpolated to U grid. | |
| 46 |
C |phiV - Varaible "phi" interpolated to V grid. | |
| 47 |
C \--------------------------------------------------------------/ |
| 48 |
C /--------------------------------------------------------------\ |
| 49 |
C |Geometry terms | |
| 50 |
C | o xa,ya,za - X, Y, Z face areas ( m.Pa, m.Pa, m**2 ). | |
| 51 |
C | o rV - 1/Volume ( m**2.Pa ) | |
| 52 |
C | o ?Grad?Op - Generic gradient operator A[XYZ][t].d/d[xyz] | |
| 53 |
C | tGradxOp: AXt/DXt, tGradyOp: AYt/DYt etc... | |
| 54 |
C \--------------------------------------------------------------/ |
| 55 |
REAL xaK (Nx+1,0:Ny) |
| 56 |
REAL yaK (0:Nx,Ny+1) |
| 57 |
REAL tGradxOp (Nx ,Ny ) |
| 58 |
REAL tGradyOp (Nx ,Ny ) |
| 59 |
REAL tGradzOpK (Nx ,Ny ) |
| 60 |
REAL tGradzOpKP1(Nx ,Ny ) |
| 61 |
REAL zaK (0:Nx,0:Ny) |
| 62 |
REAL zaKP1 (Nx,Ny) |
| 63 |
REAL rVk (Nx ,Ny ) |
| 64 |
C /--------------------------------------------------------------\ |
| 65 |
C |State variables | |
| 66 |
C | o u, v, w - X,Y,Z velocity slices. | |
| 67 |
C | o t - Potential temperature. | |
| 68 |
C | o ph - Hydrostatic pressure. | |
| 69 |
C \--------------------------------------------------------------/ |
| 70 |
REAL uK (Nx+1,0:Ny) |
| 71 |
REAL wK (0:Nx,0:Ny) |
| 72 |
REAL wkP1 (0:Nx,0:Ny) |
| 73 |
REAL vK (0:Nx,Ny+1) |
| 74 |
REAL tK (0:Nx+1,0:Ny+1) |
| 75 |
REAL phK (60) |
| 76 |
C /--------------------------------------------------------------\ |
| 77 |
C |Work space arrays | |
| 78 |
C \--------------------------------------------------------------/ |
| 79 |
REAL tmp (0:Nx+1,0:Ny+1) |
| 80 |
REAL tmp2 (0:Nx+1,0:Ny+1) |
| 81 |
C /--------------------------------------------------------------\ |
| 82 |
C |Arrays for accumulating fluxes | |
| 83 |
C | o fX, fY, fZ - Flux in X, Y and Z. | |
| 84 |
C \--------------------------------------------------------------/ |
| 85 |
REAL fX (0:Nx+1,Ny) |
| 86 |
REAL fY (Nx,0:Ny+1) |
| 87 |
REAL fZK (Nx,Ny) |
| 88 |
REAL fZKP1 (Nx,Ny) |
| 89 |
C Loop counters: |
| 90 |
INTEGER I, J |
| 91 |
C Flags: |
| 92 |
LOGICAL TOP_LAYER |
| 93 |
LOGICAL BOTTOM_LAYER |
| 94 |
C === Set up controlling flags ======================================== |
| 95 |
TOP_LAYER = K .EQ. 1 |
| 96 |
BOTTOM_LAYER = K .EQ. Nz |
| 97 |
C _D(( ' K, Nz ', K, Nz )) |
| 98 |
C _D(( ' BOTTOM_LAYER ', BOTTOM_LAYER )) |
| 99 |
C === Set up the slices =============================================== |
| 100 |
C /---------------------------------------------------------------\ |
| 101 |
C | xaK <- xa(K) | |
| 102 |
C | yaK <- ya(K) | |
| 103 |
C | za[K,KP1] <- za[(K),(K+1)] | |
| 104 |
C \---------------------------------------------------------------/ |
| 105 |
DO J=1,Ny |
| 106 |
DO I=1,Nx |
| 107 |
xaK (I,J) = xa (_I3(K,I,J)) |
| 108 |
yaK (I,J) = ya (_I3(K,I,J)) |
| 109 |
zaK (I,J) = za (_I3(K,I,J)) |
| 110 |
ENDDO |
| 111 |
ENDDO |
| 112 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
| 113 |
DO J=1,Ny |
| 114 |
DO I=1,Nx |
| 115 |
zaKP1(I,J) = za(_I3(K+1,I,J)) |
| 116 |
ENDDO |
| 117 |
ENDDO |
| 118 |
ENDIF |
| 119 |
fZK = 0. |
| 120 |
fZKP1 = 0. |
| 121 |
C /---------------------------------------------------------------\ |
| 122 |
C | tGradXOp <- 1/DX*xa | |
| 123 |
C | tGradYOp <- 1/DY*ya | |
| 124 |
C | tGradZOp[K,KP1] <- 1/DZ*ZA[K,K+1] | |
| 125 |
C \---------------------------------------------------------------/ |
| 126 |
DO J = 1, Ny |
| 127 |
DO I = 1, Nx |
| 128 |
tGradxOp (I,J) = xaK(I,J)*pDdxOp(_I3(K,I,J)) |
| 129 |
tGradyOp (I,J) = yaK(I,J)*pDdyOp(_I3(K,I,J)) |
| 130 |
tGradzOpK(I,J) = zaK(I,J)*pDdzOp(_I3(K,I,J)) |
| 131 |
ENDDO |
| 132 |
ENDDO |
| 133 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
| 134 |
DO J = 1, Ny |
| 135 |
DO I = 1, Nx |
| 136 |
tGradzOpKP1(I,J) = zaKP1(I,J)*pDdzOp(I,J,K+1) |
| 137 |
ENDDO |
| 138 |
ENDDO |
| 139 |
ENDIF |
| 140 |
C /---------------------------------------------------------------\ |
| 141 |
C | w[K,KP1] <- W[(K),(K+1)] | |
| 142 |
C \---------------------------------------------------------------/ |
| 143 |
wKP1 = 0. |
| 144 |
DO J=1,Ny |
| 145 |
DO I=1,Nx |
| 146 |
wK(I,J) = W(_I3(K,I,J)) |
| 147 |
ENDDO |
| 148 |
ENDDO |
| 149 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
| 150 |
DO J=1,Ny |
| 151 |
DO I=1,Nx |
| 152 |
wKP1(I,J) = W(_I3(K+1,I,J)) |
| 153 |
ENDDO |
| 154 |
ENDDO |
| 155 |
ENDIF |
| 156 |
C /---------------------------------------------------------------\ |
| 157 |
C | uK <- U(K) | |
| 158 |
C \---------------------------------------------------------------/ |
| 159 |
DO J=1,Ny |
| 160 |
DO I=1,Nx |
| 161 |
uK(I,J) = U(_I3(K,I,J)) |
| 162 |
ENDDO |
| 163 |
ENDDO |
| 164 |
C /---------------------------------------------------------------\ |
| 165 |
C | vK <- V(K) | |
| 166 |
C \---------------------------------------------------------------/ |
| 167 |
DO J=1,Ny |
| 168 |
DO I=1,Nx |
| 169 |
vK(I,J) = V(_I3(K,I,J)) |
| 170 |
ENDDO |
| 171 |
ENDDO |
| 172 |
C /---------------------------------------------------------------\ |
| 173 |
C | tK <- T(K) | |
| 174 |
C \---------------------------------------------------------------/ |
| 175 |
DO J=1,Ny |
| 176 |
DO I=1,Nx |
| 177 |
tK(I,J) = T(_I3(K,I,J)) |
| 178 |
ENDDO |
| 179 |
ENDDO |
| 180 |
DO J=1,Ny |
| 181 |
tK(0,J) = tK(Nx,J) |
| 182 |
tK(Nx+1,J) = tK(1 ,J) |
| 183 |
ENDDO |
| 184 |
DO I=0,Nx+1 |
| 185 |
tK(I,0) = tK(I , Ny) |
| 186 |
tK(I,Ny+1) = tK(I , 1) |
| 187 |
ENDDO |
| 188 |
C /---------------------------------------------------------------\ |
| 189 |
C | Gt Section | |
| 190 |
C |===============================================================| |
| 191 |
C | In this section Gt = GtDiffusion | |
| 192 |
C | + GtAdvection | |
| 193 |
C | + GtForcing | |
| 194 |
C \---------------------------------------------------------------/ |
| 195 |
C /---------------------------------------------------------------\ |
| 196 |
C | Initialise local terms | |
| 197 |
C | rVk <- 1/V | |
| 198 |
C \---------------------------------------------------------------/ |
| 199 |
DO J=1,Ny |
| 200 |
DO I=1,Nx |
| 201 |
rVk(I,J) = rPvol(_I3(K,I,J)) |
| 202 |
ENDDO |
| 203 |
ENDDO |
| 204 |
fZK = 0. |
| 205 |
fZKP1 = 0. |
| 206 |
C /---------------------------------------------------------------\ |
| 207 |
C | XY GtDiffusion + XY GtAdvection | |
| 208 |
C | Note: XY GtDiffusion = Laplacian + Biharmonic | |
| 209 |
C | o Laplacian ( -a2TempXY.del^2(u) ) | |
| 210 |
C | o Biharmonic ( +a4TempXY.del^4(u) ) | |
| 211 |
C \---------------------------------------------------------------/ |
| 212 |
DO J = 1, Ny |
| 213 |
DO I = 1, Nx |
| 214 |
C /-------------------------------------------------------------\ |
| 215 |
C | fX <- 1/DX*XA*Ddx{T} | |
| 216 |
C | fY <- 1/DY*YA*Ddy{T} | |
| 217 |
C \-------------------------------------------------------------/ |
| 218 |
fX(I,J) = tGradxOp(I,J)*(tK(I,J)-tK(I-1,J)) |
| 219 |
fY(I,J) = tGradyOp(I,J)*(tK(I,J)-tK(I,J-1)) |
| 220 |
ENDDO |
| 221 |
ENDDO |
| 222 |
DO J=1,Ny |
| 223 |
fX(Nx+1,J) = fX(1 ,J) |
| 224 |
ENDDO |
| 225 |
DO I=1,Nx |
| 226 |
fY(I,Ny+1) = fY(I, 1) |
| 227 |
ENDDO |
| 228 |
DO J=1,Ny |
| 229 |
DO I=1,Nx |
| 230 |
C /-------------------------------------------------------------\ |
| 231 |
C | tmp <- 1/V*(Dx{fX}+Dy{fY}) | |
| 232 |
C \-------------------------------------------------------------/ |
| 233 |
tmp(I,J)= rVk(I,J)*fX(I+1,J)-rVk(I,J)*fX(I,J) |
| 234 |
& +rVk(I,J)*fY(I,J+1)-rVk(I,J)*fY(I,J) |
| 235 |
ENDDO |
| 236 |
ENDDO |
| 237 |
tmp(0,1:Ny) = tmp(Nx,1:Ny) |
| 238 |
tmp(1:Nx,0) = tmp(1:Nx,Ny) |
| 239 |
DO J = 1,Ny |
| 240 |
DO I = 1,Nx |
| 241 |
C /-------------------------------------------------------------\ |
| 242 |
C | fX <- -k2*fX + ( Laplacian,<Line 1> ) | |
| 243 |
C | k4/DX*XA*Ddx{tmp} + ( Biharmonic,<Line 2> ) | |
| 244 |
C | XA*UBx{T} ( Advection,<Line 3> ) | |
| 245 |
C \-------------------------------------------------------------/ |
| 246 |
fX(I,J) = 0. |
| 247 |
_LPT(& -a2TempXY*fX(I,J) ) |
| 248 |
_BHT(& +a4TempXY*tGradxOp(I,J)*(tmp(I,J)-tmp(I-1,J)) ) |
| 249 |
_ADT(& +0.5*uK(I,J)*XAK(I,J)*( tK(I-1,J)+tK(I,J) ) ) |
| 250 |
C /-------------------------------------------------------------\ |
| 251 |
C | fY <- -k2*fY + ( Laplacian,<Line 1> ) | |
| 252 |
C | k4/DY*YA*Ddy{tmp} + ( Biharmonic,<Line 2> ) | |
| 253 |
C | YA*V*By{T} ( Advection,<Line 3> ) | |
| 254 |
C \-------------------------------------------------------------/ |
| 255 |
fY(I,J) = 0. |
| 256 |
_LPT(& -a2TempXY*fY(I,J) ) |
| 257 |
_BHT(& +a4TempXY*tGradyOp(I,J)*(tmp(I,J)-tmp(I,J-1)) ) |
| 258 |
_ADT(& +0.5*vK(I,J)*yaK(I,J)*( tK(I,J-1)+tK(I,J) ) ) |
| 259 |
ENDDO |
| 260 |
ENDDO |
| 261 |
C /---------------------------------------------------------------\ |
| 262 |
C | Z GtDiffusion + Z GtAdvection | |
| 263 |
C \---------------------------------------------------------------/ |
| 264 |
IF ( .NOT. TOP_LAYER ) THEN |
| 265 |
DO J = 1, Ny |
| 266 |
DO I = 1, Nx |
| 267 |
C /------------------------------------------------------------\ |
| 268 |
C | fZK <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) | |
| 269 |
C | ZA*-w*Bz{T} ( Advection,<Line 2> ) | |
| 270 |
C \------------------------------------------------------------/ |
| 271 |
fZK(I,J) = 0. |
| 272 |
_LPT(& -a2TempZ*tGradzOpK(I,J)*(T(_I3(K-1,I,J))-tK(I,J) ) ) |
| 273 |
_ADT(& +0.5*(tK(I,J)+T(_I3(K-1,I,J)))*( -wK(I,J)*ZAK(I,J) ) ) |
| 274 |
ENDDO |
| 275 |
ENDDO |
| 276 |
ELSE |
| 277 |
DO J = 1, Ny |
| 278 |
DO I = 1, Nx |
| 279 |
C /------------------------------------------------------------\ |
| 280 |
C | fZK <- ZA*-w*Bz{T} ( Advection,<Line 1> ) | |
| 281 |
C \------------------------------------------------------------/ |
| 282 |
fZK(I,J) = 0. |
| 283 |
_ADT(& +0.5*freeSurfFac*(tK(I,J)+tK(I,J))*(-wK(I,J)*ZAK(I,J) ) ) |
| 284 |
ENDDO |
| 285 |
ENDDO |
| 286 |
ENDIF |
| 287 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
| 288 |
DO J = 1, Ny |
| 289 |
DO I = 1, Nx |
| 290 |
C /------------------------------------------------------------\ |
| 291 |
C | fZKP1 <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) | |
| 292 |
C | ZA*-w*Bz{T} ( Advection,<Line 2> ) | |
| 293 |
C \------------------------------------------------------------/ |
| 294 |
fZKP1(I,J) = 0. |
| 295 |
_LPT(& -a2TempZ*tGradzOpKP1(I,J)*( tK(I,J)-T(_I3(K+1,I,J)) ) ) |
| 296 |
_ADT(& +0.5*(tK(I,J)+T(_I3(K+1,I,J)))*(-wKP1(I,J))*zaKP1(I,J) ) |
| 297 |
ENDDO |
| 298 |
ENDDO |
| 299 |
ENDIF |
| 300 |
C /---------------------------------------------------------------\ |
| 301 |
C | gT <- GtDiffusion + GtAdvection = -div(fX,fY,fZ) | |
| 302 |
C \---------------------------------------------------------------/ |
| 303 |
fX(Nx+1,1:Ny) = fX(1,1:Ny) |
| 304 |
fY(:,Ny+1) = fY(:,1) |
| 305 |
DO J = 1, Ny |
| 306 |
DO I = 1, Nx |
| 307 |
C /-------------------------------------------------------------\ |
| 308 |
C | gT <- -1/V*(Ddx{fX}+Ddy{fY}) | |
| 309 |
C \-------------------------------------------------------------/ |
| 310 |
gt(_I3(K,I,J))=-fX(I+1,J)*rVk(I,J)+fX(I,J)*rVk(I,J) |
| 311 |
& -fY(I,J+1)*rVk(I,J)+fY(I,J)*rVk(I,J) |
| 312 |
ENDDO |
| 313 |
ENDDO |
| 314 |
DO J = 1, Ny |
| 315 |
DO I = 1, Nx |
| 316 |
C /------------------------------------------------------------\ |
| 317 |
C | gT <- gT-1/Vu*fZK | |
| 318 |
C \------------------------------------------------------------/ |
| 319 |
gt(_I3(K,I,J))=gt(_I3(K,I,J))-fZK(I,J)*rVk(I,J) |
| 320 |
ENDDO |
| 321 |
ENDDO |
| 322 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
| 323 |
DO J = 1, Ny |
| 324 |
DO I = 1, Nx |
| 325 |
C /------------------------------------------------------------\ |
| 326 |
C | gT <- gT+1/Vu*fZP1 | |
| 327 |
C \------------------------------------------------------------/ |
| 328 |
gt(_I3(K,I,J))=gt(_I3(K,I,J))+fZKP1(I,J)*rVk(I,J) |
| 329 |
ENDDO |
| 330 |
ENDDO |
| 331 |
ENDIF |
| 332 |
C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MAXVAL(Gt )', MAXVAL(Gt))) |
| 333 |
C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MINVAL(Gt )', MINVAL(Gt))) |
| 334 |
C /---------------------------------------------------------------\ |
| 335 |
C | gT <- gT + GtForcing | |
| 336 |
C \---------------------------------------------------------------/ |
| 337 |
IF ( TOP_LAYER ) THEN |
| 338 |
DO J =1,Ny |
| 339 |
DO I=1, Nx |
| 340 |
Cdbg gT(_I3(K,I,J)) = gT(_I3(K,I,J)) + HEAT(I,J) |
| 341 |
ENDDO |
| 342 |
ENDDO |
| 343 |
ENDIF |
| 344 |
|
| 345 |
C Final time derivative from AB2. |
| 346 |
gt(_I3(K,1:Nx,1:Ny)) = Gt(_I3(K,1:Nx,1:Ny))*pMask(_I3(K,:,:)) |
| 347 |
tmp = 0. |
| 348 |
tmp2 = 0. |
| 349 |
tmp (1:Nx,1:Ny) = gtNM1(_I3(K,:,:)) |
| 350 |
tmp2(1:Nx,1:Ny) = gt (_I3(K,:,:)) |
| 351 |
|
| 352 |
C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(GtNM1)', MAXVAL(GtNm1) )) |
| 353 |
C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(Gt)', MAXVAL(Gt) )) |
| 354 |
|
| 355 |
gtNM1(_I3(K,:,:)) = gt(_I3(K,:,:)) |
| 356 |
tmp2=(1.5+abEps)*tmp2-(0.5+abEps)*tmp |
| 357 |
gt(_I3(K,1:Nx,1:Ny))=tmp2(1:Nx,1:Ny) |
| 358 |
C |
| 359 |
RETURN |
| 360 |
END |
Parent Directory
| 
Revision Log
| 
Revision Graph