*IF DEF,SEAICE ORH1F305.445
C ******************************COPYRIGHT****************************** GTS2F400.4267
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.4268
C GTS2F400.4269
C Use, duplication or disclosure of this code is subject to the GTS2F400.4270
C restrictions as set forth in the contract. GTS2F400.4271
C GTS2F400.4272
C Meteorological Office GTS2F400.4273
C London Road GTS2F400.4274
C BRACKNELL GTS2F400.4275
C Berkshire UK GTS2F400.4276
C RG12 2SZ GTS2F400.4277
C GTS2F400.4278
C If no contract has been raised with this copy of the code, the use, GTS2F400.4279
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.4280
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.4281
C Modelling at the above address. GTS2F400.4282
C ******************************COPYRIGHT****************************** GTS2F400.4283
C GTS2F400.4284
C*LL ICEADVEC.4
CLL SUBROUTINE ICE_ADVECT ICEADVEC.5
CLL ------------------- ICEADVEC.6
CLL ICEADVEC.7
CLL DYNAMIC SEA ICE MODEL SUBROUTINE TO ADVECT ICE THICKNESS, ICEADVEC.8
CLL COMPACTNESS AND SNOW DEPTH ON ARAKAWA C GRID USING UPSTREAM ICEADVEC.9
CLL DIFFERENCING. ICEADVEC.10
CLL ICEADVEC.11
CLL IT CAN BE COMPILED BY CFT77, BUT DOES NOT CONFORM TO ANSI ICEADVEC.12
CLL FORTRAN77 STANDARDS, BECAUSE OF THE INLINE COMMENTS. ICEADVEC.13
CLL ICEADVEC.14
CLL ALL QUANTITIES IN THIS ROUTINE ARE IN S.I. UNITS UNLESS ICEADVEC.15
CLL OTHERWISE STATED. ICEADVEC.16
CLL ICEADVEC.17
CLL 15/09/93 Advection scheme now conserves ice mass taking full ICEADVEC.18
CLL account of spherical geometry. JFT ICEADVEC.19
CLL Modified 08/10/93 T.C.Johns ICEADVEC.20
CLL Converted to nupdate format for compatibility with vn3.2. ICEADVEC.21
CLL ICEADVEC.22
CLL WRITTEN BY J.F.THOMSON (18/02/93) ICEADVEC.23
CLL ICEADVEC.24
CLL MODEL MODIFICATION HISTORY FROM INSERTION INTO UM 3.3: ICEADVEC.25
CLL VERSION DATE ICEADVEC.26
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.4750
! 4.3 01.02.97 Include further SWAPBOUNDS for mpp. R.Hill ORH3F403.60
CLL ICEADVEC.27
CLL THIS ROUTINE FORMS PART OF SYSTEM COMPONENT P4. ICEADVEC.28
CLL ADHERES TO THE STANDARDS OF DOCUMENTATION PAPER 4, VERSION 1. ICEADVEC.29
CLLEND--------------------------------------------------------------- ICEADVEC.30
C*L ICEADVEC.31
subroutine ice_advect( 1,7ICEADVEC.32
*CALL ARGOINDX 
ORH7F402.237
& uice ICEADVEC.33
&,vice ICEADVEC.34
&,aice ICEADVEC.35
&,hice ICEADVEC.36
&,hsnow ICEADVEC.37
&,hmask ICEADVEC.38
&,umask ICEADVEC.39
&,icy ICEADVEC.40
&,imt,imtm2,jmt,jmtm1 ORH3F405.98
&,phi ICEADVEC.42
&,phit ICEADVEC.43
&,dphi ICEADVEC.44
&,dlambda ICEADVEC.45
&,radius ICEADVEC.46
&,dtts ICEADVEC.47
& ) ICEADVEC.48
C ICEADVEC.49
implicit none ICEADVEC.50
C ICEADVEC.51
*CALL OTIMER ORH1F305.4751
*CALL TYPOINDX ORH7F402.238
integer ICEADVEC.52
& imt ! number of tracer columns. ICEADVEC.53
&,imtm2 ! number of tracer columns - 2 ORH6F402.109
&,jmt ! number of tracer rows. ICEADVEC.54
&,jmtm1 ! number of velocity rows. ICEADVEC.55
real ICEADVEC.56
& uice(imt,jmtm1) ! in zonal sea ice velocity. ICEADVEC.57
&,vice(imt,jmtm1) ! in meridional sea ice velocity. ICEADVEC.58
&,phi(jmt) ! in latitude of uv rows in radians. ICEADVEC.59
&,phit(jmt) ! in latitude of mass rows in radians. ICEADVEC.60
&,aice(imt,jmt) ! inout ice compactness. ICEADVEC.61
&,hice(imt,jmt) ! inout ice thickness. ICEADVEC.62
&,hsnow(imt,jmt) ! inout snow depth. ICEADVEC.63
real ICEADVEC.64
& radius ! in radius of the earth in metres. ICEADVEC.65
&,dphi(jmt) ! in meridional grid spacing in radians. ICEADVEC.66
&,dlambda(imt) ! in zonal grid spacing in radians. ICEADVEC.67
&,dtts ! in timestep in seconds. ICEADVEC.68
real ICEADVEC.69
& hmask(imt,jmt) ! in 1.0 for ha land 0.0 for sea. ICEADVEC.70
&,umask(imt,jmtm1) ! in 1.0 for uv land 0.0 for sea. ICEADVEC.71
logical ICEADVEC.72
& icy(imt,jmt) ! in true if ice is present. ICEADVEC.73
C ICEADVEC.74
C variables local to this subroutine are now defined ICEADVEC.75
C ICEADVEC.76
real ICEADVEC.77
& hice_init(imt,jmt) ! initial ice thickness. ICEADVEC.78
&,aice_init(imt,jmt) ! initial ice compactness. ICEADVEC.79
&,hsno_init(imt,jmt) ! initial snow depth. ICEADVEC.80
&,uwork(imt,jmtm1),vwork(imt,jmtm1) ICEADVEC.81
real ICEADVEC.82
& um,vm,ump,vmp ICEADVEC.83
&,qinx,qinxp,qiny,qinyp ICEADVEC.84
&,qax,qaxp,qay,qayp ICEADVEC.85
&,qsx,qsxp,qsy,qsyp ICEADVEC.86
&,qijh,qija,qijs ICEADVEC.87
&,area,facey,faceyp,facex,facexp ICEADVEC.88
C ICEADVEC.89
logical ICEADVEC.90
& lcheck ICEADVEC.91
C ICEADVEC.92
integer ICEADVEC.93
& i ORH6F402.111
&,j ICEADVEC.98
C* ICEADVEC.99
! start executable code ORH9F402.207
ORH9F402.208
if (L_OTIMER) call timer('iceadvec',3) ORH1F305.4752
C ICEADVEC.108
*IF DEF,MPP ORH4F402.135
C===================================================================== ORH4F402.136
C CALL TO SWAPBOUNDS FOR HALO UPDATE IN MPP VERSION ORH4F402.137
C===================================================================== ORH4F402.138
ORH4F402.139
CALL SWAPBOUNDS(UICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ORH3F403.61
CALL SWAPBOUNDS(VICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ORH3F403.62
ORH3F403.63
ORH4F402.143
*ENDIF ORH4F402.144
ORH4F402.145
C Copy initial ice thickness, compactness and snow depth to workspace ICEADVEC.109
C ICEADVEC.110
do j=J_1,J_jmt ORH9F402.209
do i=1,imt ICEADVEC.112
hice_init(i,j) = hice(i,j) ICEADVEC.113
aice_init(i,j) = aice(i,j) ICEADVEC.114
hsno_init(i,j) = hsnow(i,j) ICEADVEC.115
end do ICEADVEC.116
end do ICEADVEC.117
*IF DEF,MPP ORH4F402.100
C===================================================================== ORH4F402.101
C CALL TO SWAPBOUNDS FOR HALO UPDATE IN MPP VERSION ORH4F402.102
C===================================================================== ORH4F402.103
ORH4F402.104
CALL SWAPBOUNDS(HICE_INIT,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH4F402.105
ORH4F402.106
CALL SWAPBOUNDS(AICE_INIT,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH4F402.107
ORH4F402.108
CALL SWAPBOUNDS(HSNO_INIT,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ORH4F402.109
ORH4F402.110
*ENDIF ORH4F402.111
ORH4F402.112
C ICEADVEC.118
C Loop over grid advecting primary variables. ICEADVEC.119
C ICEADVEC.120
ORH9F402.210
! Changed the following loop from j=1,jmtm2 ORH9F402.211
! to j=2,jmtm1 to avoid referencing variables outside ORH9F402.212
! halo regions in MPP code. Hence all refs to ORH9F402.213
! j+1 become j, j+2 becomes j+1 and j becomes j-1 ORH9F402.214
do j=J_2,J_jmtm1 ORH9F402.215
do i=1,imtm2 ICEADVEC.122
lcheck = (icy(i,j).or.icy(i+1,j).or.icy(i+1,j-1) ORH9F402.216
& .or.icy(i+1,j+1).or.icy(i+2,j)) ORH9F402.217
if ( hmask(i+1,j).gt.0.1 .and. lcheck ) then ORH9F402.218
C ICEADVEC.126
area =radius**2*dlambda(i+1)*(sin(phi(j))-sin(phi(j-1))) ORH9F402.219
facex=radius*dphi(j) ORH9F402.220
facexp = facex ICEADVEC.129
facey = radius*cos(phi(j-1))*dlambda(i+1) ORH9F402.221
faceyp = radius*cos(phi(j))*dlambda(i+1) ORH9F402.222
C ICEADVEC.132
um = uice(i,j-1) ORH9F402.223
vm = vice(i,j-1) ORH9F402.224
ump = uice(i+1,j-1) ORH9F402.225
vmp = vice(i,j) ORH9F402.226
C ICEADVEC.137
if (um.ge.0.) qinx = um*hice_init(i,j)*hmask(i,j) ORH9F402.227
if (um.lt.0.) qinx = um*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.228
if (vm.ge.0.) qiny = vm*hice_init(i+1,j-1)*hmask(i+1,j-1) ORH9F402.229
if (vm.lt.0.) qiny = vm*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.230
if (ump.ge.0.) qinxp=-ump*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.231
if (ump.lt.0.) qinxp=-ump*hice_init(i+2,j)*hmask(i+2,j) ORH9F402.232
if (vmp.ge.0.) qinyp=-vmp*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.233
if (vmp.lt.0.) qinyp=-vmp*hice_init(i+1,j+1)*hmask(i+1,j+1) ORH9F402.234
C ICEADVEC.146
if (um.ge.0.) qax = um*aice_init(i,j)*hmask(i,j) ORH9F402.235
if (um.lt.0.) qax = um*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.236
if (vm.ge.0.) qay = vm*aice_init(i+1,j-1)*hmask(i+1,j-1) ORH9F402.237
if (vm.lt.0.) qay = vm*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.238
if (ump.ge.0.) qaxp =-ump*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.239
if (ump.lt.0.) qaxp =-ump*aice_init(i+2,j)*hmask(i+2,j) ORH9F402.240
if (vmp.ge.0.) qayp =-vmp*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.241
if (vmp.lt.0.) qayp =-vmp*aice_init(i+1,j+1)*hmask(i+1,j+1) ORH9F402.242
C ICEADVEC.155
if (um.ge.0.) qsx = qax*hsno_init(i,j) ORH9F402.243
if (um.lt.0.) qsx = qax*hsno_init(i+1,j) ORH9F402.244
if (vm.ge.0.) qsy = qay*hsno_init(i+1,j-1) ORH9F402.245
if (vm.lt.0.) qsy = qay*hsno_init(i+1,j) ORH9F402.246
if (ump.ge.0.) qsxp = qaxp*hsno_init(i+1,j) ORH9F402.247
if (ump.lt.0.) qsxp = qaxp*hsno_init(i+2,j) ORH9F402.248
if (vmp.ge.0.) qsyp = qayp*hsno_init(i+1,j) ORH9F402.249
if (vmp.lt.0.) qsyp = qayp*hsno_init(i+1,j+1) ORH9F402.250
C ICEADVEC.164
qijh = hice(i+1,j)*hmask(i+1,j) ORH9F402.251
qija = aice(i+1,j)*hmask(i+1,j) ORH9F402.252
qijs = hsnow(i+1,j)*qija ORH9F402.253
C ICEADVEC.168
hice(i+1,j) = qijh + ( qinx*facex + qinxp*facexp + qiny*facey ORH9F402.254
& + qinyp*faceyp ) * dtts/area ICEADVEC.170
aice(i+1,j) = qija + ( qax*facex + qaxp*facexp + qay*facey ORH9F402.255
& + qayp*faceyp ) * dtts/area ICEADVEC.172
if (aice(i+1,j).gt.0.0) ORH9F402.256
& hsnow(i+1,j)= ( qijs + ( qsx*facex + qsxp*facexp + qsy*facey ORH9F402.257
& + qsyp*faceyp ) * dtts/area )/aice(i+1,j) ORH9F402.258
C ICEADVEC.176
endif ICEADVEC.177
end do ICEADVEC.178
end do ICEADVEC.179
! ORH9F402.259
*IF DEF,MPP ORH9F402.260
IF (JFIN.GE.JMT_GLOBAL.AND.JST.LE.JMT_GLOBAL) THEN ORH9F402.261
*ENDIF ORH9F402.262
ORH9F402.263
! Special code for row JMT_GLOBAL ORH9F402.264
! ORH9F402.265
j = J_jmt ORH9F402.266
do i=1,imtm2 ORH9F402.267
lcheck = (icy(i,j).or.icy(i+1,j).or.icy(i+1,j-1) ORH9F402.268
& .or.icy(i+2,j)) ORH9F402.269
if (hmask(i+1,j).gt.0.1 .and. lcheck) then ORH9F402.270
C ICEADVEC.180
area=radius**2*dlambda(i+1)*(sin(phi(j))-sin(phi(j-1))) ORH9F402.271
facex = radius*dphi(j) ORH9F402.272
facexp = facex ICEADVEC.191
facey = radius*cos(phi(j-1))*dlambda(i+1) ORH9F402.273
C ICEADVEC.193
um = uice(i,j-1) ORH9F402.274
vm = vice(i,j-1) ORH9F402.275
ump = uice(i+1,j-1) ORH9F402.276
vmp = 0. ICEADVEC.197
C ICEADVEC.198
if (um.ge.0.) qinx = um*hice_init(i,j)*hmask(i,j) ORH9F402.277
if (um.lt.0.) qinx = um*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.278
if (vm.ge.0.) qiny = vm*hice_init(i+1,j-1)*hmask(i+1,j-1) ORH9F402.279
if (vm.lt.0.) qiny = vm*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.280
if (ump.ge.0.) qinxp =-ump*hice_init(i+1,j)*hmask(i+1,j) ORH9F402.281
if (ump.lt.0.) qinxp =-ump*hice_init(i+2,j)*hmask(i+2,j) ORH9F402.282
qinyp = 0. ICEADVEC.205
C ICEADVEC.206
if (um.ge.0.) qax = um*aice_init(i,j)*hmask(i,j) ORH9F402.283
if (um.lt.0.) qax = um*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.284
if (vm.ge.0.) qay = vm*aice_init(i+1,j-1)*hmask(i+1,j-1) ORH9F402.285
if (vm.lt.0.) qay = vm*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.286
if (ump.ge.0.) qaxp =-ump*aice_init(i+1,j)*hmask(i+1,j) ORH9F402.287
if (ump.lt.0.) qaxp =-ump*aice_init(i+2,j)*hmask(i+2,j) ORH9F402.288
qayp = 0.0 ICEADVEC.213
C ICEADVEC.214
if (um.ge.0.) qsx = qax*hsno_init(i,j) ORH9F402.289
if (um.lt.0.) qsx = qax*hsno_init(i+1,j) ORH9F402.290
if (vm.ge.0.) qsy = qay*hsno_init(i+1,j-1) ORH9F402.291
if (vm.lt.0.) qsy = qay*hsno_init(i+1,j) ORH9F402.292
if (ump.ge.0.) qsxp = qaxp*hsno_init(i+1,j) ORH9F402.293
if (ump.lt.0.) qsxp = qaxp*hsno_init(i+2,j) ORH9F402.294
qsyp = 0. ICEADVEC.221
C ICEADVEC.222
qijh = hice(i+1,j)*hmask(i+1,j) ORH9F402.295
qija = aice(i+1,j)*hmask(i+1,j) ORH9F402.296
qijs = hsnow(i+1,j)*qija ORH9F402.297
C ICEADVEC.226
hice(i+1,j) = qijh + ( qinx*facex + qinxp*facexp + qiny*facey ) ORH9F402.298
& * dtts/area ICEADVEC.228
aice(i+1,j) = qija + ( qax*facex + qaxp*facexp + qay*facey ) ORH9F402.299
& * dtts/area ICEADVEC.230
hsnow(i+1,j)= ( qijs + ( qsx*facex + qsxp*facexp + qsy*facey ) ORH9F402.300
& * dtts/area )/aice(i+1,j) ORH9F402.301
endif ICEADVEC.233
end do ICEADVEC.234
*IF DEF,MPP ORH9F402.302
ENDIF ! If this PE deals with row JMT_GLOBAL ORH9F402.303
*ENDIF ORH9F402.304
ORH9F402.305
C ICEADVEC.237
if (L_OTIMER) call timer('iceadvec',4) ORH1F305.4753
return ICEADVEC.241
end ICEADVEC.242
*ENDIF ICEADVEC.243