subroutine slab_temp_advect( 1,4SLBTADV.25
& L1,u_field,landmask SJC1F400.297
&,icols,jrows,jrowsm1,lglobal SJC1F400.298
&,delta_lat,delta_long,timestep,a,dz1 SJC1F400.299
&,usea,vsea,tmask,opensea SLBTADV.28
&,cos_p_latitude,sec_p_latitude SJC1F400.300
&,cos_u_latitude,sin_u_latitude SJC1F400.301
&,slabtemp,wtsfc,wtbase SJC1F400.302
& ) SLBTADV.31
C SLBTADV.32
implicit none SLBTADV.33
C SLBTADV.34
integer SLBTADV.35
& L1 ! in points in p field SJC1F400.303
&,u_field ! in points in u field SJC1F400.304
&,icols ! in number of columns E-W SJC1F400.305
&,jrows ! in number of rows N-S SLBTADV.37
&,jrowsm1 ! in number of rows N-S - 1 SLBTADV.38
logical SLBTADV.39
& lglobal ! in true for a global model SLBTADV.40
&,landmask(icols,jrows) ! in land-sea mask (p-grid) SJC1F400.306
&,opensea(icols,jrows) ! in true if box is open sea SLBTADV.41
real SLBTADV.42
& delta_lat ! in EW grid spacing in degrees SLBTADV.43
&,delta_long ! in NS grid spacing in degrees SLBTADV.44
&,timestep ! in timestep in seconds SLBTADV.45
&,a ! in radius of earth in metres SLBTADV.46
&,dz1 ! in slab ocean thickness (m) SJC1F400.307
&,usea(icols,jrowsm1) ! in zonal surface current (M/S) SLBTADV.47
&,vsea(icols,jrowsm1) ! in meri surface current (M/S) SLBTADV.48
&,tmask(icols,jrows) !in mask 1.0 for opensea 0.0 la/ice SLBTADV.49
real SLBTADV.50
& cos_p_latitude(icols,jrows) ! in cos(latitude) on p grid SLBTADV.51
&,cos_u_latitude(icols,jrowsm1) ! in cos(latitude) on uv grid SLBTADV.52
&,sin_u_latitude(icols,jrowsm1) ! in sin(latitude) on uv grid SLBTADV.53
&,sec_p_latitude(icols,jrowsm1) ! in sec(latitude) on p grid SJC1F400.308
&,slabtemp(icols,jrows) ! inout temperatue of slab ocean C SLBTADV.54
&,wtsfc(icols,jrows) ! inout w * surface slab temp SJC1F400.309
&,wtbase(icols,jrows) ! inout w * base slab temp SJC1F400.310
C SLBTADV.55
C Variables local to this subroutine are now defined SJC1F400.311
C SJC1F400.312
EXTERNAL ZONM SJC1F400.313
C SLBTADV.57
real SLBTADV.58
& slabt_init(icols,jrows) ! initial slab temperature SLBTADV.59
&,wbase(icols,jrows) ! vertical velocity at base SJC1F400.314
&,wsfc ! vertical velocity at surface =0 SJC1F400.315
&,divuv(icols,jrows) ! divergence of u and v SJC1F400.316
&,um ! usea at i,j SLBTADV.60
&,ump ! usea at i+1,j SLBTADV.61
&,vm ! vsea at i,j SLBTADV.62
&,vmp ! vsea at i,j+i SLBTADV.63
&,wm ! wsea at surface SJC1F400.317
&,wmp ! wsea at base SJC1F400.318
&,tinx ! advection across facex SLBTADV.64
&,tinxp ! advection across facexp SLBTADV.65
&,tiny ! advection across facey SLBTADV.66
&,tinyp ! advection across faceyp SLBTADV.67
&,tinz ! advection across face top SJC1F400.319
&,tinzp ! advection across face bottom SJC1F400.320
&,tinzsum ! sum of vertical advection comps SJC1F400.321
&,tinzwt ! total weight of upwelling points SJC1F400.322
&,tiny_pole ! advection across facey at pole SLBTADV.68
&,tinyp_pole ! advection across faceyp at pole SLBTADV.69
&,area ! grid box area SLBTADV.70
&,fractarea ! area*0.125 SLBTADV.71
&,facex ! length of west side of box SLBTADV.72
&,facexp ! length of east side of box SLBTADV.73
&,facey ! length of north side of box SLBTADV.74
&,faceyp ! length of south side of box SLBTADV.75
&,wtbpwts ! wtbase + wtsfc for conservation SJC1F400.323
&,p_levels ! no of p levels = 1 SJC1F400.324
&,latitude_step_inverse ! 1 / latitude increment SJC1F400.325
&,longitude_step_inverse ! 1 / longitude increment SJC1F400.326
&,smask(icols,jrows) ! Sea mask (p-grid) for zonal mean SJC1F400.327
&,s_pmass(icols,jrows) ! dummy wgt for surf(p_grid) SJC1F400.328
&,z_slabtemp(jrows) ! zonal mean slab temp SJC1F400.329
&,dtarea ! timestep/area SJC1F400.330
&,dtdz ! timestep/slab depth SJC1F400.331
&,cosplat ! cos_p_latitude SJC1F400.332
&,coslimit ! cos(latitude limit for vertadv) SJC1F400.333
parameter ( coslimit = 0.642 ) ! limit is 50 degrees SJC1F400.334
C SLBTADV.76
integer SLBTADV.77
& icolsm1 ! icols - 1 SLBTADV.78
&,icolsm2 ! icols - 2 SLBTADV.79
&,jrowsm2 ! jrows - 2 SLBTADV.80
&,i,j,i1,i2,ii ! loop counts SJC1F400.335
&,spts(jrows) ! No of sea points/row SJC1F400.336
C SJC1F400.337
LOGICAL SJC1F400.338
& lspts(jrows) ! Marks rows with sea pts SJC1F400.339
C* SLBTADV.82
C start executable code SLBTADV.83
C SLBTADV.84
*IF DEF,TIMER SLBTADV.85
call timer
('slbtadv',3) SLBTADV.86
*ENDIF SLBTADV.87
icolsm1 = icols-1 SLBTADV.88
icolsm2 = icols-2 SLBTADV.89
jrowsm2 = jrows-2 SLBTADV.90
tinyp_pole = 0.0 SLBTADV.91
tiny_pole = 0.0 SLBTADV.92
p_levels = 1 SJC1F400.340
tinzsum = 0.0 SJC1F400.341
tinzwt = 0.0 SJC1F400.342
dtdz = timestep / dz1 SJC1F400.343
latitude_step_inverse = 1. / delta_lat SJC1F400.344
longitude_step_inverse = 1. / delta_long SJC1F400.345
C SLBTADV.95
C 1. Calculate zonal mean temperature SJC1F400.346
C SJC1F400.347
C 1.1 Set up masks for weighted sums SJC1F400.348
C SJC1F400.349
DO j=1,jrows SJC1F400.350
DO i=1,icols SJC1F400.351
IF (.not. LANDMASK(I,j)) THEN SJC1F400.352
SMASK(I,j) = 1.0 SJC1F400.353
else SJC1F400.354
SMASK(I,j) = 0.0 SJC1F400.355
ENDIF SJC1F400.356
END DO SJC1F400.357
END DO SJC1F400.358
C SJC1F400.359
C 1.2 Calculate no of sea points on row-by-row basis SJC1F400.360
C and set logical array to denote active sea rows SJC1F400.361
C SJC1F400.362
DO j=1,jROWS SJC1F400.363
SPTS(j) = 0 SJC1F400.364
DO i=1,icols SJC1F400.365
SPTS(j) = SPTS(j) + SMASK(I,j) SJC1F400.366
end do SJC1F400.367
if (spts(j) .gt.0) then SJC1F400.368
LSPTS(j) = .true. SJC1F400.369
else SJC1F400.370
LSPTS(j) = .false. SJC1F400.371
endif SJC1F400.372
end do SJC1F400.373
C SJC1F400.374
C 1.3 Set dummy weighting for surface variables to one SJC1F400.375
C SJC1F400.376
DO j=1,jrows SJC1F400.377
DO i=1,icols SJC1F400.378
S_PMASS(i,j) = 1.0 SJC1F400.379
end do SJC1F400.380
end do SJC1F400.381
C SJC1F400.382
C 1.4 Mass weighted zonal mean slabtemp on p grid SJC1F400.383
C SJC1F400.384
call zonm(slabtemp,z_slabtemp SJC1F400.385
&,smask,s_pmass,lspts,icols,jrows) SJC1F400.386
C SJC1F400.387
C 2. Calculate vertical velocity SJC1F400.388
C SJC1F400.389
C with boundary condition wmp=0 at surface SJC1F400.390
wsfc = 0.0 SJC1F400.391
C SJC1F400.392
C 2.1 Divergence of horizontal velocity SJC1F400.393
C SJC1F400.394
call div_calc(usea,vsea,u_field,L1,p_levels, SJC1F400.395
& icols,sec_p_latitude,cos_u_latitude, SJC1F400.396
& latitude_step_inverse,longitude_step_inverse,1,divuv) SJC1F400.397
C SJC1F400.398
C 2.2 Vertical velocity = -dz1 * div (u) SJC1F400.399
C wbase positive downwards SJC1F400.400
C SJC1F400.401
do j=1,jrowsm2 SJC1F400.402
do i=1,icols SJC1F400.403
wbase(i,j+1) = -1. * dz1 * divuv(i,j+1) SJC1F400.404
end do SJC1F400.405
end do SJC1F400.406
C SJC1F400.407
C 3. Copy initial slab temperature to workspace SJC1F400.408
C SLBTADV.97
do j=1,jrows SLBTADV.98
do i=1,icols SLBTADV.99
slabt_init(i,j) = slabtemp(i,j) SLBTADV.100
end do SLBTADV.101
end do SLBTADV.102
C SLBTADV.103
C 4. Conservation of vertical advection SJC1F400.409
C SJC1F400.410
C Make vertical advection of slab temperature conservative. SJC1F400.411
C As wtbase and wtsfc are calculated add to get total sum, SJC1F400.412
C to conserve this sum needs to be zero. SJC1F400.413
C Total made zero by adding difference onto upwelling points. SJC1F400.414
C This is done using 5 iterations. SJC1F400.415
C Vertical advection only applied between latitude limit (50 N/S) SJC1F400.416
C SJC1F400.417
C 4.1 Calculate wtbase and wtsfc and total vertical advection SJC1F400.418
C SLBTADV.105
do j=1,jrowsm2 SLBTADV.106
cosplat = cos_p_latitude(1,j+1) SJC1F400.419
C SLBTADV.112
do i=1,icols SLBTADV.113
i1 = i+1 SLBTADV.114
i2 = i+2 SLBTADV.115
if (i1.gt.icols) i1 = i1-icols SLBTADV.116
if (i2.gt.icols) i2 = i2-icols SLBTADV.117
C SLBTADV.118
if ( tmask(i1,j+1) .gt. 0.1 ) then SJC1F400.420
C SJC1F400.421
wm = wsfc SJC1F400.422
wmp = wbase(i1,j+1) SJC1F400.423
C SJC1F400.424
if ( cosplat .ge. coslimit ) then SJC1F400.425
C SJC1F400.426
C Downwelling advection terms SJC1F400.427
if (wmp .gt. 0.) then SJC1F400.428
wtsfc(i1,j+1) = -wmp*slabt_init(i1,j+1) SJC1F400.429
wtbase(i1,j+1) = 0.0 SJC1F400.430
tinzsum = tinzsum +wtsfc(i1,j+1)*cosplat SJC1F400.431
endif SJC1F400.432
C SJC1F400.433
C No vertical advection at these points SJC1F400.434
if (wmp .eq. 0.) then SJC1F400.435
wtbase(i1,j+1) = 0.0 SJC1F400.436
wtsfc(i1,j+1) = 0.0 SJC1F400.437
endif SJC1F400.438
C SJC1F400.439
C Upwelling advection terms SJC1F400.440
if (wmp.lt.0.) then SJC1F400.441
wtbase(i1,j+1) = -wmp*z_slabtemp(j+1) SJC1F400.442
wtsfc(i1,j+1) = 0.0 SJC1F400.443
tinzwt = tinzwt + 1.0*cosplat SJC1F400.444
tinzsum = tinzsum + wtbase(i1,j+1)*cosplat SJC1F400.445
endif SJC1F400.446
C SJC1F400.447
else SJC1F400.448
C These points are polewards of area of vertical advection. SJC1F400.449
wtsfc(i1,j+1) = 0.0 SJC1F400.450
wtbase(i1,j+1) = 0.0 SJC1F400.451
endif SJC1F400.452
endif SJC1F400.453
end do SJC1F400.454
end do SJC1F400.455
C SJC1F400.456
C 4.2 tinsum = weighted sum of vertical advection terms SJC1F400.457
C want to make this zero SJC1F400.458
C so add it to each upwelling point SJC1F400.459
C and divide by weights of upwelling points SJC1F400.460
C SJC1F400.461
wtbpwts = tinzsum / tinzwt SJC1F400.462
C SJC1F400.463
C 4.3 Make tinzsum zero using 5 iterations SJC1F400.464
C SJC1F400.465
do ii=1,5 SJC1F400.466
tinzwt = 0.0 SJC1F400.467
tinzsum = 0.0 SJC1F400.468
cosplat = 0.0 SJC1F400.469
C SJC1F400.470
do j=1,jrowsm2 SJC1F400.471
cosplat = cos_p_latitude(1,j+1) SJC1F400.472
C SJC1F400.473
do i=1,icols SJC1F400.474
i1 = i+1 SJC1F400.475
i2 = i+2 SJC1F400.476
if (i1.gt.icols) i1 = i1-icols SJC1F400.477
if (i2.gt.icols) i2 = i2-icols SJC1F400.478
C SJC1F400.479
if ( tmask(i1,j+1) .gt. 0.1 ) then SJC1F400.480
wmp = wbase(i1,j+1) SJC1F400.481
C SJC1F400.482
if ( cosplat .ge. coslimit ) then SJC1F400.483
C SJC1F400.484
C Change upwelling points by wtbpwts SJC1F400.485
C Add upwelling to total and calculate total weight SJC1F400.486
if (wmp.lt.0.) then SJC1F400.487
wtbase(i1,j+1) = wtbase(i1,j+1) - wtbpwts SJC1F400.488
tinzwt = tinzwt + 1.0 * cosplat SJC1F400.489
tinzsum = tinzsum + wtbase(i1,j+1)*cosplat SJC1F400.490
endif SJC1F400.491
C SJC1F400.492
C Add downwelling to total SJC1F400.493
if (wmp.gt.0.) then SJC1F400.494
tinzsum = tinzsum + wtsfc(i1,j+1)*cosplat SJC1F400.495
endif SJC1F400.496
C SJC1F400.497
endif SJC1F400.498
endif SJC1F400.499
enddo SJC1F400.500
enddo SJC1F400.501
C SJC1F400.502
C Recalculate difference to add to upweeling points SJC1F400.503
wtbpwts = tinzsum / tinzwt SJC1F400.504
enddo SJC1F400.505
C SJC1F400.506
C 5. Loop over grid advecting slab temperature SJC1F400.507
C SJC1F400.508
cosplat =0.0 SJC1F400.509
facex = ABS(a * delta_lat) SJC1F400.510
facexp = facex SJC1F400.511
C SJC1F400.512
C 5.1 All but polar rows SJC1F400.513
C SJC1F400.514
do j=1,jrowsm2 SJC1F400.515
area = ABS(a**2 * delta_long * (sin_u_latitude(1,j+1) SJC1F400.516
& -sin_u_latitude(1,j))) SJC1F400.517
fractarea = area * 0.125 SJC1F400.518
dtarea = timestep / area SJC1F400.519
facey = ABS(a*cos_u_latitude(1,j)*delta_long) SJC1F400.520
faceyp = ABS(a*cos_u_latitude(1,j+1)*delta_long) SJC1F400.521
cosplat = cos_p_latitude(1,j+1) SJC1F400.522
C SJC1F400.523
do i=1,icols SJC1F400.524
i1 = i+1 SJC1F400.525
i2 = i+2 SJC1F400.526
if (i1.gt.icols) i1 = i1-icols SJC1F400.527
if (i2.gt.icols) i2 = i2-icols SJC1F400.528
C SJC1F400.529
if ( tmask(i1,j+1) .gt. 0.1 ) then SJC1F400.530
C SLBTADV.120
um = usea(i,j) SLBTADV.121
vm = vsea(i,j) SLBTADV.122
ump = usea(i1,j) SLBTADV.123
vmp = vsea(i,j+1) SLBTADV.124
wm = wsfc SJC1F400.531
wmp = wbase(i1,j+1) SJC1F400.532
C SLBTADV.125
C Vertical components SJC1F400.533
if (wmp .gt. 0.) then SJC1F400.534
tinzp = wtsfc(i1,j+1) SJC1F400.535
endif SJC1F400.536
if (wmp .eq. 0.) then SJC1F400.537
tinzp = 0.0 SJC1F400.538
endif SJC1F400.539
if (wmp .lt. 0.) then SJC1F400.540
tinzp = wtbase(i1,j+1) SJC1F400.541
endif SJC1F400.542
tinz = wm SJC1F400.543
C SJC1F400.544
C Horizontal components SJC1F400.545
if (um.ge.0.) tinx = um*slabt_init(i,j+1)*tmask(i,j+1) SLBTADV.126
if (um.lt.0.) tinx = um*slabt_init(i1,j+1)*tmask(i1,j+1) SLBTADV.127
if (vm.ge.0.) tiny =-vm*slabt_init(i1,j+1)*tmask(i1,j+1) SLBTADV.128
if (vm.lt.0.) tiny =-vm*slabt_init(i1,j)*tmask(i1,j) SLBTADV.129
if (ump.ge.0.) tinxp =-ump*slabt_init(i1,j+1)*tmask(i1,j+1) SLBTADV.130
if (ump.lt.0.) tinxp =-ump*slabt_init(i2,j+1)*tmask(i2,j+1) SLBTADV.131
if (vmp.ge.0.) tinyp = vmp*slabt_init(i1,j+2)*tmask(i1,j+2) SLBTADV.132
if (vmp.lt.0.) tinyp = vmp*slabt_init(i1,j+1)*tmask(i1,j+1) SLBTADV.133
C SLBTADV.134
C Polar rows SJC1F400.546
if (j.eq.1) tiny_pole = tiny_pole SLBTADV.135
& -tiny*facey*timestep/fractarea SLBTADV.136
if (j.eq.jrowsm2) tinyp_pole = tinyp_pole SLBTADV.137
& -tinyp*faceyp*timestep/fractarea SLBTADV.138
C SLBTADV.139
C Advection equation SJC1F400.547
slabtemp(i1,j+1) = slabtemp(i1,j+1) SLBTADV.140
& + ( tinx*facex + tinxp*facexp SLBTADV.141
& + tiny*facey + tinyp*faceyp ) * dtarea SJC1F400.548
& + ( tinz + tinzp ) * dtdz SJC1F400.549
C SLBTADV.143
endif SLBTADV.144
end do SLBTADV.145
end do SLBTADV.146
C SLBTADV.147
C 5.2 Special code for polar rows. SJC1F400.550
C SLBTADV.149
tiny_pole = tiny_pole/icols SLBTADV.150
tinyp_pole = tinyp_pole/icols SLBTADV.151
C SLBTADV.152
C First row SJC1F400.551
do i=1,icols SLBTADV.154
i1 = i+1 SLBTADV.155
i2 = i+2 SLBTADV.156
if (i1.gt.icols) i1 = i1-icols SLBTADV.157
if (i2.gt.icols) i2 = i2-icols SLBTADV.158
if (tmask(i1,1).gt.0.1) then SLBTADV.159
slabtemp(i1,1) = slabtemp(i1,1) + tiny_pole SLBTADV.160
endif SLBTADV.161
end do SLBTADV.162
C Last row SJC1F400.552
do i=1,icols SLBTADV.164
i1 = i+1 SLBTADV.165
i2 = i+2 SLBTADV.166
if (i1.gt.icols) i1 = i1-icols SLBTADV.167
if (i2.gt.icols) i2 = i2-icols SLBTADV.168
if (tmask(i1,jrows).gt.0.1) then SLBTADV.169
slabtemp(i1,jrows) = slabtemp(i1,jrows) + tinyp_pole SLBTADV.170
endif SLBTADV.171
end do SLBTADV.172
C SLBTADV.173
*IF DEF,TIMER SLBTADV.174
call timer('slbtadv',4) SLBTADV.175
*ENDIF SLBTADV.176
return SLBTADV.177
end SLBTADV.178
*ENDIF SLBTADV.179