*IF DEF,OCEAN CALCRL1A.2
C *****************************COPYRIGHT****************************** CALCRL1A.3
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. CALCRL1A.4
C CALCRL1A.5
C Use, duplication or disclosure of this code is subject to the CALCRL1A.6
C restrictions as set forth in the contract. CALCRL1A.7
C CALCRL1A.8
C Meteorological Office CALCRL1A.9
C London Road CALCRL1A.10
C BRACKNELL CALCRL1A.11
C Berkshire UK CALCRL1A.12
C RG12 2SZ CALCRL1A.13
C CALCRL1A.14
C If no contract has been raised with this copy of the code, the use, CALCRL1A.15
C duplication or disclosure of it is strictly prohibited. Permission CALCRL1A.16
C to do so must first be obtained in writing from the Head of Numerical CALCRL1A.17
C Modelling at the above address. CALCRL1A.18
C ******************************COPYRIGHT****************************** CALCRL1A.19
!+ Routine to calculate rigid lid surface pressure CALCRL1A.20
! Subroutine Interface: CALCRL1A.21
CALCRL1A.22
SUBROUTINE CALC_RLIDP( 2,12CALCRL1A.23
*CALL ARGSIZE
CALCRL1A.24
*CALL ARGOCALL
CALCRL1A.25
*CALL ARGOINDX
CALCRL1A.26
& ICODE,CMESSAGE CALCRL1A.27
&,ITT,ZU,ZV,PTD CALCRL1A.28
&,rlsrfp) CALCRL1A.29
CALCRL1A.30
IMPLICIT NONE CALCRL1A.31
CALCRL1A.32
! Purpose: CALCRL1A.33
! To calculate the rigid lid pressure over the model domain CALCRL1A.34
! CALCRL1A.35
! Method: CALCRL1A.36
! Calculates sea surface pressure gradient from vertically averaged CALCRL1A.37
! forcing and streamfunction tendency. CALCRL1A.38
! Integrates sea surface pressure gradient over the domain and CALCRL1A.39
! subtracts area average to obtain the surface pressure to within a CALCRL1A.40
! constant. Divides by 1000 to yield the sea surface height. CALCRL1A.41
! CALCRL1A.42
! Current Code Owner: R.Forbes CALCRL1A.43
! CALCRL1A.44
! History: CALCRL1A.45
! Version Date Comment CALCRL1A.46
! ------- ---- ------- CALCRL1A.47
! 4.4 01/10/97 New code R.Forbes CALCRL1A.48
!LL 4.5 17/09/98 Update calls to timer, required because of GPB8F405.80
!LL new barrier inside timer. P.Burton GPB8F405.81
! CALCRL1A.49
! Code Description: CALCRL1A.50
! Language: FORTRAN 77 + common extensions. CALCRL1A.51
! This code is written to UMDP3 v6 programming standards. CALCRL1A.52
! CALCRL1A.53
! ---------------------------------------------------------------------- CALCRL1A.54
! Global variables: CALCRL1A.55
*CALL UMSCALAR
CALCRL1A.56
CALCRL1A.57
! Subroutine arguments: CALCRL1A.58
CALCRL1A.59
*CALL OARRYSIZ
CALCRL1A.60
*CALL TYPSIZE
CALCRL1A.61
*CALL TYPOINDX
PXORDER.8
*CALL TYPOCALL
CALCRL1A.62
CALCRL1A.64
! Scalar arguments with intent(in): CALCRL1A.65
INTEGER ITT ! time step CALCRL1A.66
CALCRL1A.67
! Array arguments with intent(in): CALCRL1A.68
REAL ZU(IMT,JMT) ! U-component of vertically averaged forcing CALCRL1A.69
REAL ZV(IMT,JMT) ! V-component of vertically averaged forcing CALCRL1A.70
REAL PTD(IMT_STREAM, JMT_STREAM) ! streamfunction tendency CALCRL1A.71
CALCRL1A.72
! Array arguments with intent(out): CALCRL1A.73
REAL rlsrfp(IMT,JMT) ! sea surface height (cm) on TS-grid CALCRL1A.74
CALCRL1A.75
! ErrorStatus: CALCRL1A.76
INTEGER ICODE ! Return code : 0 Normal Exit : >0 Error CALCRL1A.77
CHARACTER*(80) CMESSAGE ! Error message if return code >0 CALCRL1A.78
CALCRL1A.79
! Local Parameters: CALCRL1A.80
*CALL CNTLOCN
CALCRL1A.81
*CALL OTIMER
CALCRL1A.82
*CALL C_MDI
CALCRL1A.83
*CALL GCCOM
CALCRL1A.84
CALCRL1A.85
! Local Scalars: CALCRL1A.86
INTEGER i,j ! loop indices CALCRL1A.87
INTEGER ndry ! number of land points CALCRL1A.88
INTEGER nwet ! number of ocean points CALCRL1A.89
INTEGER info CALCRL1A.90
REAL surfpx,surfpy ! rl surface pressure gradient CALCRL1A.91
REAL dubdt,dvbdt ! } CALCRL1A.92
REAL diago3,diago4 ! } variables for rigid lid surface CALCRL1A.93
REAL fxa,f2,f3 ! } pressure gradient calculation CALCRL1A.94
REAL atosp,detmr ! } CALCRL1A.95
REAL r2dtuv ! reciprocal of baroclinic timestep CALCRL1A.96
CALCRL1A.97
! Local Arrays: CALCRL1A.98
REAL spx(imt,jmt) ! x-gradient in surface pressure on uv-grid CALCRL1A.99
REAL spy(imt,jmt) ! y-gradient in surface pressure on uv-grid CALCRL1A.100
LOGICAL wet(imt,jmt) ! T=> ocean point on ts-grid CALCRL1A.101
LOGICAL wetu(imt,jmt) ! T=> ocean point on uv-grid CALCRL1A.102
CALCRL1A.103
! Local Arrays on mpp "global" grid CALCRL1A.104
REAL spx_gl(imt,jmt_global) CALCRL1A.105
REAL spy_gl(imt,jmt_global) CALCRL1A.106
REAL dyu_gl(jmt_global) ! y grid spacing on uv-grid CALCRL1A.107
REAL dyt_gl(jmt_global) ! y grid spacing on ts-grid CALCRL1A.108
REAL cs_gl(jmt_global) ! cos(lat) on uv-grid CALCRL1A.109
REAL cst_gl(jmt_global) ! cos(lat) on ts-grid CALCRL1A.110
REAL rlsrfp_gl(imt,jmt_global) ! rigid-lid surface pressure CALCRL1A.111
LOGICAL wet_gl(imt,jmt_global) ! T=> ocean point on ts-grid CALCRL1A.112
LOGICAL wetu_gl(imt,jmt_global) ! T=> ocean point on uv-grid CALCRL1A.113
CALCRL1A.114
! Function & Subroutine calls: CALCRL1A.115
EXTERNAL RLP_LINEINT CALCRL1A.116
CALCRL1A.117
!- End of header CALCRL1A.118
! ---------------------------------------------------------------------- CALCRL1A.119
CALCRL1A.120
IF (L_OTIMER) CALL TIMER
('CALC_RLP',103) GPB8F405.82
CALCRL1A.122
! Initialize arrays CALCRL1A.123
DO j=j_1,j_jmt CALCRL1A.124
DO i=1,imt CALCRL1A.125
spx(i,j) = 0.0 CALCRL1A.126
spy(i,j) = 0.0 CALCRL1A.127
rlsrfp(i,j) = 0.0 CALCRL1A.128
wet(i,j) = .false. CALCRL1A.129
wetu(i,j) = .false. CALCRL1A.130
ENDDO CALCRL1A.131
ENDDO CALCRL1A.132
CALCRL1A.133
! On mixing time steps "ptd" has been multiplied by a factor of 2 CALCRL1A.134
! and the time step has to be adjusted also. CALCRL1A.135
CALCRL1A.136
fxa = 1.0 CALCRL1A.137
IF (MIX.EQ.1) THEN CALCRL1A.138
fxa = 0.5 CALCRL1A.139
ENDIF CALCRL1A.140
CALCRL1A.141
! --------------------------------------------------------------------- CALCRL1A.142
! Compute surface pressure gradient on uv-grid CALCRL1A.143
! Add in the external mode part of d/dt for the momentum balance CALCRL1A.144
! Also add in external mode part of the implicit coriolis term CALCRL1A.145
! Also set mask "wet" on TS-grid CALCRL1A.146
! --------------------------------------------------------------------- CALCRL1A.147
CALCRL1A.148
r2dtuv = 1.0/c2dtuv CALCRL1A.149
DO j=j_1,j_jmtm1 CALCRL1A.150
CALCRL1A.151
atosp = acor*2.0*omega*sine(j) CALCRL1A.152
f2 = atosp*c2dtuv CALCRL1A.153
detmr = 1.0/(1.0+f2*f2) CALCRL1A.154
IF (.NOT.(L_ONOCLIN)) f3 = c2dtuv/c2dtsf CALCRL1A.155
CALCRL1A.156
DO i=1,imtm1 CALCRL1A.157
CALCRL1A.158
wet(i,j) = fkmp(i,j).GT.0 CALCRL1A.159
wetu(i,j) = fkmq(i,j).GT.0 CALCRL1A.160
CALCRL1A.161
IF (wetu(i,j)) THEN CALCRL1A.162
! This is an ocean grid point CALCRL1A.163
CALCRL1A.164
IF (.NOT.(L_ONOCLIN))THEN CALCRL1A.165
! Section 31: Barotropic solution CALCRL1A.166
diago3 = fxa*(ptd(i+1,j+1)-ptd(i ,j)) CALCRL1A.167
diago4 = fxa*(ptd(i ,j+1)-ptd(i+1,j)) CALCRL1A.168
dubdt = (diago3+diago4)*dyu2r(j)*hr(i,j) CALCRL1A.169
dvbdt = (diago3-diago4)*dxu2r(i)*hr(i,j)*csr(j) CALCRL1A.170
surfpx = r2dtuv*( dubdt + f3*zu(i,j) + f2*dvbdt) CALCRL1A.171
surfpy = r2dtuv*(-dvbdt + f3*zv(i,j) + f2*dubdt) CALCRL1A.172
ELSE CALCRL1A.173
! Section 30: No barotropic solution CALCRL1A.174
surfpx = r2dtuv*detmr*(zu(i,j) + zv(i,j)*f2) CALCRL1A.175
surfpy = r2dtuv*detmr*(zv(i,j) - zu(i,j)*f2) CALCRL1A.176
ENDIF ! L_ONOCLIN CALCRL1A.177
CALCRL1A.178
spx(i,j) = surfpx CALCRL1A.179
spy(i,j) = surfpy CALCRL1A.180
CALCRL1A.181
ENDIF ! wetu CALCRL1A.182
CALCRL1A.183
ENDDO CALCRL1A.184
ENDDO CALCRL1A.185
CALCRL1A.186
! Apply cyclic boundary condition CALCRL1A.187
CALCRL1A.188
IF (L_OCYCLIC) THEN CALCRL1A.189
DO j=j_1,j_jmtm1 CALCRL1A.190
spx(1,j) = spx(imtm1,j) CALCRL1A.191
spx(imt,j) = spx(2,j) CALCRL1A.192
spy(1,j) = spy(imtm1,j) CALCRL1A.193
spy(imt,j) = spy(2,j) CALCRL1A.194
wet(1,j) = wet(imtm1,j) CALCRL1A.195
wet(imt,j) = wet(2,j) CALCRL1A.196
wetu(1,j) = wetu(imtm1,j) CALCRL1A.197
wetu(imt,j)= wetu(2,j) CALCRL1A.198
ENDDO CALCRL1A.199
ENDIF CALCRL1A.200
CALCRL1A.201
! Find number of wet and dry TS-points CALCRL1A.202
CALCRL1A.203
ndry = 0 CALCRL1A.204
nwet = 0 CALCRL1A.205
DO j=j_1,j_jmt CALCRL1A.206
DO i=1,imt CALCRL1A.207
IF (wet(i,j)) THEN CALCRL1A.208
nwet = nwet + 1 CALCRL1A.209
ELSE CALCRL1A.210
ndry = ndry + 1 CALCRL1A.211
ENDIF CALCRL1A.212
ENDDO CALCRL1A.213
ENDDO CALCRL1A.214
CALCRL1A.215
! --------------------------------------------------------------------- CALCRL1A.216
! Perform integration to obtain surface pressure field CALCRL1A.217
! --------------------------------------------------------------------- CALCRL1A.218
CALCRL1A.219
*IF DEF,MPP CALCRL1A.220
! Gather local fields to form global field CALCRL1A.221
CALCRL1A.222
CALL O_SMARTPASS
(1,1,dyu(J_1),dyu_gl CALCRL1A.223
& ,jfin-jst+1,jmt_global,jst,2) CALCRL1A.224
CALCRL1A.225
CALL O_SMARTPASS
(1,1,dyt(J_1),dyt_gl CALCRL1A.226
& ,jfin-jst+1,jmt_global,jst,2) CALCRL1A.227
CALCRL1A.228
CALL O_SMARTPASS
(1,1,cs(J_1),cs_gl CALCRL1A.229
& ,jfin-jst+1,jmt_global,jst,2) CALCRL1A.230
CALCRL1A.231
CALL O_SMARTPASS
(1,1,cst(J_1),cst_gl CALCRL1A.232
& ,jfin-jst+1,jmt_global,jst,2) CALCRL1A.233
CALCRL1A.234
CALL GATHER_FIELD
(spx,spx_gl, CALCRL1A.235
& imt,jmt,imt,jmt_global, CALCRL1A.236
& 0,GC_ALLGROUP,info) CALCRL1A.237
CALCRL1A.238
CALL GATHER_FIELD
(spy,spy_gl, CALCRL1A.239
& imt,jmt,imt,jmt_global, CALCRL1A.240
& 0,GC_ALLGROUP,info) CALCRL1A.241
CALCRL1A.242
CALL GATHER_FIELD
(wet,wet_gl, CALCRL1A.243
& imt,jmt,imt,jmt_global, CALCRL1A.244
& 0,GC_ALLGROUP,info) CALCRL1A.245
CALCRL1A.246
CALL GATHER_FIELD
(wetu,wetu_gl, CALCRL1A.247
& imt,jmt,imt,jmt_global, CALCRL1A.248
& 0,GC_ALLGROUP,info) CALCRL1A.249
CALCRL1A.250
CALL GC_ISUM(
1,O_NPROC,info,nwet) CALCRL1A.251
CALL GC_ISUM(
1,O_NPROC,info,ndry) CALCRL1A.252
CALCRL1A.253
! Switch to single processor CALCRL1A.254
IF (o_mype .EQ. 0) THEN CALCRL1A.255
CALCRL1A.256
*ENDIF CALCRL1A.257
CALCRL1A.258
! Perform line integral CALCRL1A.259
CALL RLP_LINEINT
(icode,cmessage,imt,jmt_global,ndry,nwet, CALCRL1A.260
& dxu,dyu_gl,dxt,dyt_gl,cs_gl,cst_gl, CALCRL1A.261
& spx_gl,spy_gl,rlsrfp_gl,l_ocyclic,wet_gl,wetu_gl ) CALCRL1A.262
CALCRL1A.263
*IF DEF,MPP CALCRL1A.264
ENDIF ! am I PE 0 ? CALCRL1A.265
CALCRL1A.266
! Scatter global field to local fields CALCRL1A.267
CALL SCATTER_FIELD
(rlsrfp,rlsrfp_gl, CALCRL1A.268
& imt,jmt,imt,jmt_global, CALCRL1A.269
& 0,GC_ALLGROUP,info) CALCRL1A.270
CALCRL1A.271
*ENDIF CALCRL1A.272
CALCRL1A.273
999 CONTINUE CALCRL1A.274
CALCRL1A.275
IF (L_OTIMER) CALL TIMER
('CALC_RLP',104) GPB8F405.83
CALCRL1A.277
RETURN CALCRL1A.278
END CALCRL1A.279
! --------------------------------------------------------------------- CALCRL1A.280
!+ Routine to perform line integral of a 2-D gradient field CALCRL1A.281
! Subroutine Interface: CALCRL1A.282
CALCRL1A.283
SUBROUTINE RLP_LINEINT(icode,cmessage, 1CALCRL1A.284
& imt,jmt,ndry,nwet,dxu,dyu,dxt,dyt,cs,cst, CALCRL1A.285
& spx,spy,rlsrfp,l_ocyclic,wet,wetu) CALCRL1A.286
CALCRL1A.287
IMPLICIT NONE CALCRL1A.288
CALCRL1A.289
! Purpose: CALCRL1A.290
! perform line integral of a 2-D gradient field CALCRL1A.291
! CALCRL1A.292
! Method: CALCRL1A.293
! CALCRL1A.294
! Current Code Owner: R.Forbes CALCRL1A.295
! CALCRL1A.296
! History: CALCRL1A.297
! Version Date Comment CALCRL1A.298
! ------- ---- ------- CALCRL1A.299
! 4.4 01/10/97 New code R.Forbes CALCRL1A.300
! CALCRL1A.301
! Code Description: CALCRL1A.302
! Language: FORTRAN 77 + common extensions. CALCRL1A.303
! This code is written to UMDP3 v6 programming standards. CALCRL1A.304
! CALCRL1A.305
! ---------------------------------------------------------------------- CALCRL1A.306
! Array arguments with intent(in): CALCRL1A.307
INTEGER imt,jmt ! grid dimensions CALCRL1A.308
INTEGER ndry,nwet ! no. of land and ocean points on ts-grid CALCRL1A.309
REAL dxu(imt) ! x grid spacing on uv-grid CALCRL1A.310
REAL dyu(jmt) ! y grid spacing on uv-grid CALCRL1A.311
REAL dxt(imt) ! x grid spacing on ts-grid CALCRL1A.312
REAL dyt(jmt) ! y grid spacing on ts-grid CALCRL1A.313
REAL cs(jmt) ! cos(lat) on uv-grid CALCRL1A.314
REAL cst(jmt) ! cos(lat) on ts-grid CALCRL1A.315
REAL spx(imt,jmt) ! x-gradient in surface pressure on uv-grid CALCRL1A.316
REAL spy(imt,jmt) ! y-gradient in surface pressure on uv-grid CALCRL1A.317
LOGICAL wet(imt,jmt) ! T=> ocean point on ts-grid CALCRL1A.318
LOGICAL wetu(imt,jmt) ! T=> ocean point on uv-grid CALCRL1A.319
LOGICAL L_OCYCLIC ! T=> cyclic grid CALCRL1A.320
CALCRL1A.321
! Array arguments with intent(out): CALCRL1A.322
REAL rlsrfp(IMT,JMT) ! sea surface height (cm) on TS-grid CALCRL1A.323
CALCRL1A.324
! ErrorStatus: CALCRL1A.325
INTEGER ICODE ! Return code : 0 Normal Exit : >0 Error CALCRL1A.326
CHARACTER*(80) CMESSAGE ! Error message if return code >0 CALCRL1A.327
CALCRL1A.328
! Local Scalars: CALCRL1A.329
INTEGER maxpass ! maximum no. of passes for fill-in loop CALCRL1A.330
INTEGER npass ! counter for number of passes CALCRL1A.331
INTEGER gridpass ! loop index over chequerboard grids CALCRL1A.332
INTEGER i,j,ii,jj ! loop indices CALCRL1A.333
INTEGER nfilled ! number of filled points CALCRL1A.334
INTEGER noldfilled ! number of filled points CALCRL1A.335
INTEGER ilon0,jlat0 ! } CALCRL1A.336
INTEGER istrti,jstrti ! } start and end grid points CALCRL1A.337
INTEGER istopi,jstopi ! } for line integral CALCRL1A.338
INTEGER jrev ! } CALCRL1A.339
REAL pavg ! area mean surface pressure CALCRL1A.340
REAL wetarea ! total area covered by ocean CALCRL1A.341
CALCRL1A.342
! Local Arrays: CALCRL1A.343
LOGICAL filled(imt,jmt) CALCRL1A.344
CALCRL1A.345
! --------------------------------------------------------------------- CALCRL1A.346
! Initialize parameters CALCRL1A.347
maxpass = 200 CALCRL1A.348
CALCRL1A.349
! Initialize arrays CALCRL1A.350
DO j=1,jmt CALCRL1A.351
DO i=1,imt CALCRL1A.352
rlsrfp(i,j) = 0.0 CALCRL1A.353
filled(i,j) = .false. CALCRL1A.354
ENDDO CALCRL1A.355
ENDDO CALCRL1A.356
CALCRL1A.357
! --------------------------------------------------------------------- CALCRL1A.358
! make first pass over domain integrating until land is reached CALCRL1A.359
! (integration on TS-grid) CALCRL1A.360
! --------------------------------------------------------------------- CALCRL1A.361
CALCRL1A.362
! Hardwire start point for surface pressure integration CALCRL1A.363
! for Tropical Pacific grid CALCRL1A.364
! Note: these parameters should be set in the user-interface CALCRL1A.365
! or calculated from the land sea mask CALCRL1A.366
ilon0 = 80 CALCRL1A.367
jlat0 = 60 CALCRL1A.368
CALCRL1A.369
! Loop over two grids for integration CALCRL1A.370
DO gridpass = 1,2 CALCRL1A.371
CALCRL1A.372
! make sure starting out in ocean CALCRL1A.373
IF (.NOT.wet(ilon0,jlat0)) THEN CALCRL1A.374
WRITE(6,9901) ilon0,jlat0 CALCRL1A.375
9901 FORMAT(' ERROR: START POINT (',i3,',',i3,') NOT IN OCEAN') CALCRL1A.376
ICODE=1 CALCRL1A.377
CMESSAGE='ERROR: START POINT NOT IN OCEAN' CALCRL1A.378
GOTO 999 CALCRL1A.379
ENDIF CALCRL1A.380
CALCRL1A.381
! start integration with pressure equal to zero (arbitrary) CALCRL1A.382
rlsrfp(ilon0,jlat0) = 0.0 CALCRL1A.383
filled(ilon0,jlat0) = .true. CALCRL1A.384
c CALCRL1A.385
C --------------------------------------------------------------------- CALCRL1A.386
C integrate north-east starting from ilon0,jlat0 CALCRL1A.387
C --------------------------------------------------------------------- CALCRL1A.388
jstopi = jlat0 CALCRL1A.389
i = ilon0 CALCRL1A.390
IF (jlat0+1 .LE. jmt-1) THEN CALCRL1A.391
DO j=jlat0+1,jmt-1 CALCRL1A.392
i = i + 1 CALCRL1A.393
IF ((i .LE. imt-1) .AND. wet(i,j) .AND. wetu(i-1,j-1)) THEN CALCRL1A.394
istopi = i CALCRL1A.395
jstopi = j CALCRL1A.396
rlsrfp(i,j) = rlsrfp(i-1,j-1) + CALCRL1A.397
& spx(i-1,j-1)*dxu(i-1)*cs(j-1) + spy(i-1,j-1)*dyu(j-1) CALCRL1A.398
filled(i,j) = .true. CALCRL1A.399
ENDIF CALCRL1A.400
ENDDO CALCRL1A.401
ENDIF CALCRL1A.402
CALCRL1A.403
! --------------------------------------------------------------------- CALCRL1A.404
! integrate south-west starting from ilon0,jlat0 CALCRL1A.405
! --------------------------------------------------------------------- CALCRL1A.406
jstrti = jlat0 CALCRL1A.407
i=ilon0 CALCRL1A.408
IF (jlat0-1 .ge. 2) THEN CALCRL1A.409
DO jrev=2,jlat0-1 CALCRL1A.410
j = jlat0 - jrev + 1 CALCRL1A.411
i = i - 1 CALCRL1A.412
IF ((i .GE. 2) .AND. wet(i,j) .AND. wetu(i+1,j+1)) THEN CALCRL1A.413
istrti = i CALCRL1A.414
jstrti = j CALCRL1A.415
rlsrfp(i,j) = rlsrfp(i+1,j+1) - CALCRL1A.416
& spx(i,j)*dxu(i)*cs(j) - spy(i,j)*dyu(j) CALCRL1A.417
filled(i,j) = .true. CALCRL1A.418
ENDIF CALCRL1A.419
ENDDO CALCRL1A.420
ENDIF CALCRL1A.421
CALCRL1A.422
! --------------------------------------------------------------------- CALCRL1A.423
! integrate south-east from initial long. until land is reached CALCRL1A.424
! --------------------------------------------------------------------- CALCRL1A.425
ii = istrti - 1 CALCRL1A.426
DO jj=jstrti,jstopi CALCRL1A.427
CALCRL1A.428
ii = ii + 1 CALCRL1A.429
i = ii CALCRL1A.430
DO j=jj-1,2,-1 CALCRL1A.431
i = i + 1 CALCRL1A.432
IF ((i .LE. imt-1) .AND. wet(i,j) .AND. wetu(i-1,j+1)) THEN CALCRL1A.433
rlsrfp(i,j) = rlsrfp(i-1,j+1) + CALCRL1A.434
& spx(i-1,j)*dxu(i-1)*cs(j) - spy(i-1,j)*dyu(j) CALCRL1A.435
filled(i,j) = .true. CALCRL1A.436
ENDIF CALCRL1A.437
ENDDO CALCRL1A.438
CALCRL1A.439
! --------------------------------------------------------------------- CALCRL1A.440
! integrate north-west from initial long. until land is reached CALCRL1A.441
! --------------------------------------------------------------------- CALCRL1A.442
i = ii CALCRL1A.443
DO j=jj+1,jmt-1 CALCRL1A.444
i = i - 1 CALCRL1A.445
IF ((i .GE. 2) .AND. wet(i,j) .AND. wetu(i+1,j-1)) THEN CALCRL1A.446
rlsrfp(i,j) = rlsrfp(i+1,j-1) - CALCRL1A.447
& spx(i,j-1)*dxu(i)*cs(j-1) + spy(i,j-1)*dyu(j-1) CALCRL1A.448
filled(i,j) = .true. CALCRL1A.449
ENDIF CALCRL1A.450
ENDDO CALCRL1A.451
CALCRL1A.452
ENDDO CALCRL1A.453
CALCRL1A.454
CALCRL1A.455
IF (L_OCYCLIC) THEN CALCRL1A.456
DO j=1,jmt-1 CALCRL1A.457
rlsrfp(1,j) = rlsrfp(imt-1,j) CALCRL1A.458
rlsrfp(imt,j) = rlsrfp(2,j) CALCRL1A.459
filled(1,j) = filled(imt-1,j) CALCRL1A.460
filled(imt,j) = filled(2,j) CALCRL1A.461
ENDDO CALCRL1A.462
ENDIF CALCRL1A.463
CALCRL1A.464
nfilled = 0 CALCRL1A.465
DO j=1,jmt CALCRL1A.466
DO i=1,imt CALCRL1A.467
IF (wet(i,j) .AND. filled(i,j)) CALCRL1A.468
& nfilled = nfilled + 1 CALCRL1A.469
ENDDO CALCRL1A.470
ENDDO CALCRL1A.471
CALCRL1A.472
! --------------------------------------------------------------------- CALCRL1A.473
! pass over entire domain looking for unfilled wet points CALCRL1A.474
! --------------------------------------------------------------------- CALCRL1A.475
npass = 0 CALCRL1A.476
noldfilled = 0 CALCRL1A.477
CALCRL1A.478
DO WHILE (nfilled .LT. nwet .AND. npass .LT. maxpass .AND. CALCRL1A.479
& nfilled .GT. noldfilled) CALCRL1A.480
CALCRL1A.481
noldfilled = nfilled CALCRL1A.482
CALCRL1A.483
DO j=1,jmt CALCRL1A.484
DO i=1,imt CALCRL1A.485
IF (wet(i,j) .AND. (.NOT. filled(i,j)) ) THEN CALCRL1A.486
! this point needs to be filled, look at surrounding points CALCRL1A.487
CALCRL1A.488
! look south-east CALCRL1A.489
IF ((i .ne. imt) .AND. (j .ne. 1) .AND. CALCRL1A.490
& wet(i+1,j-1) .AND. filled(i+1,j-1) .AND. CALCRL1A.491
& wetu(i,j-1)) THEN CALCRL1A.492
rlsrfp(i,j) = rlsrfp(i+1,j-1) - CALCRL1A.493
& spx(i,j-1)*dxu(i)*cs(j-1) + spy(i,j-1)*dyu(j-1) CALCRL1A.494
filled(i,j) = .true. CALCRL1A.495
CALCRL1A.496
! look south-west CALCRL1A.497
ELSEIF ((i .ne. 1) .AND. (j .ne. 1) .AND. CALCRL1A.498
& wet(i-1,j-1) .AND. filled(i-1,j-1) .AND. CALCRL1A.499
& wetu(i-1,j-1)) THEN CALCRL1A.500
rlsrfp(i,j) = rlsrfp(i-1,j-1) + CALCRL1A.501
& spx(i-1,j-1)*dxu(i-1)*cs(j-1) + spy(i-1,j-1)*dyu(j-1) CALCRL1A.502
filled(i,j) = .true. CALCRL1A.503
CALCRL1A.504
! look north-west CALCRL1A.505
ELSEIF ((i .ne. 1) .AND. (j .ne. jmt) .AND. CALCRL1A.506
& wet(i-1,j+1) .AND. filled(i-1,j+1) .AND. CALCRL1A.507
& wetu(i-1,j)) THEN CALCRL1A.508
rlsrfp(i,j) = rlsrfp(i-1,j+1) + CALCRL1A.509
& spx(i-1,j)*dxu(i-1)*cs(j) - spy(i-1,j)*dyu(j) CALCRL1A.510
filled(i,j) = .true. CALCRL1A.511
CALCRL1A.512
! look north-east CALCRL1A.513
ELSEIF ((i .ne. imt) .AND. (j .ne. jmt) .AND. CALCRL1A.514
& wet(i+1,j+1) .AND. filled(i+1,j+1) .AND. CALCRL1A.515
& wetu(i,j)) THEN CALCRL1A.516
rlsrfp(i,j) = rlsrfp(i+1,j+1) - CALCRL1A.517
& spx(i,j)*dxu(i)*cs(j) - spy(i,j)*dyu(j) CALCRL1A.518
filled(i,j) = .true. CALCRL1A.519
CALCRL1A.520
ENDIF CALCRL1A.521
ENDIF CALCRL1A.522
ENDDO CALCRL1A.523
ENDDO CALCRL1A.524
CALCRL1A.525
npass = npass + 1 CALCRL1A.526
CALCRL1A.527
IF (L_OCYCLIC) THEN CALCRL1A.528
DO j=1,jmt-1 CALCRL1A.529
rlsrfp(1,j) = rlsrfp(imt-1,j) CALCRL1A.530
rlsrfp(imt,j) = rlsrfp(2,j) CALCRL1A.531
filled(1,j) = filled(imt-1,j) CALCRL1A.532
filled(imt,j) = filled(2,j) CALCRL1A.533
ENDDO CALCRL1A.534
ENDIF CALCRL1A.535
CALCRL1A.536
nfilled = 0 CALCRL1A.537
DO j=1,jmt CALCRL1A.538
DO i=1,imt CALCRL1A.539
IF (wet(i,j) .AND. filled(i,j)) nfilled = nfilled + 1 CALCRL1A.540
ENDDO CALCRL1A.541
ENDDO CALCRL1A.542
CALCRL1A.543
ENDDO ! on while CALCRL1A.544
CALCRL1A.545
CALCRL1A.546
! if not all points filled after maxpass, write out indices CALCRL1A.547
IF ((nwet/2)-1 .GT. nfilled) THEN CALCRL1A.548
DO j=1,jmt CALCRL1A.549
DO i=1,imt CALCRL1A.550
IF (wet(i,j) .AND. (.NOT. filled(i,j))) CALCRL1A.551
& WRITE(6,9905) i,j CALCRL1A.552
9905 FORMAT(' UNFILLED POINT AT I,J = ',2i5) CALCRL1A.553
ENDDO CALCRL1A.554
ENDDO CALCRL1A.555
ICODE=1 CALCRL1A.556
CMESSAGE='RLIDSRFP: ERROR - UNFILLED POINTS IN INTEGRATION' CALCRL1A.557
GOTO 999 CALCRL1A.558
ENDIF CALCRL1A.559
CALCRL1A.560
CALCRL1A.561
! compute area mean pressure CALCRL1A.562
pavg = 0.0 CALCRL1A.563
wetarea = 0.0 CALCRL1A.564
DO j=1,jmt CALCRL1A.565
DO i=1,imt CALCRL1A.566
IF (wet(i,j) .AND. filled(i,j)) THEN CALCRL1A.567
pavg = pavg + rlsrfp(i,j)*dxt(i)*cst(j)*dyt(j) CALCRL1A.568
wetarea =wetarea + dxt(i)*cst(j)*dyt(j) CALCRL1A.569
ENDIF CALCRL1A.570
ENDDO CALCRL1A.571
ENDDO CALCRL1A.572
pavg = pavg/wetarea CALCRL1A.573
CALCRL1A.574
CALCRL1A.575
! normalize pressure CALCRL1A.576
DO j=1,jmt CALCRL1A.577
DO i=1,imt CALCRL1A.578
IF (wet(i,j) .AND. filled(i,j)) CALCRL1A.579
& rlsrfp(i,j) = (rlsrfp(i,j) - pavg) CALCRL1A.580
ENDDO CALCRL1A.581
ENDDO CALCRL1A.582
ilon0 = ilon0 + 1 CALCRL1A.583
CALCRL1A.584
ENDDO ! on gridpass CALCRL1A.585
CALCRL1A.586
999 CONTINUE CALCRL1A.587
CALCRL1A.588
RETURN CALCRL1A.589
END CALCRL1A.590
*ENDIF CALCRL1A.591