SUBROUTINE icefreedr( 1,14ICEFREED.36
*CALL ARGOINDX 
ICEFREED.37
& imt,imtm1,imtm2,jmt,jmtm1,jmtm2, ICEFREED.38
& rhoice,rhowater,rhosnow, ICEFREED.39
& cw,hicestop,hiceslow, ICEFREED.40
& coriolis, ICEFREED.41
& icy,ocean,aice,hice,hsnow, ICEFREED.42
& ucurrent,vcurrent, ICEFREED.43
& wsx_ice,wsy_ice, ICEFREED.44
& uice,vice, ICEFREED.45
& isx,isy, ICEFREED.46
& twaterx,twatery, ICEFREED.47
& mfx,mfy, ICEFREED.48
& inisx,inisy, ICEFREED.49
& radius,dphi,dlambda,phit ICEFREED.50
&) ICEFREED.51
IMPLICIT none ICEFREED.52
*CALL CNTLOCN ICEFREED.53
*CALL OARRYSIZ ICEFREED.54
*CALL TYPOINDX ICEFREED.55
! Declare arguments ICEFREED.56
INTEGER ICEFREED.57
& imt ! IN # columns ICEFREED.58
&,imtm1 ICEFREED.59
&,imtm2 ICEFREED.60
&,jmt ! IN # tracer rows ICEFREED.61
&,jmtm1 ! IN # velocity rows ICEFREED.62
&,jmtm2 ICEFREED.63
REAL ICEFREED.64
& rhoice, ! IN SI density of seaice ICEFREED.65
& rhowater, ! IN SI density of seawater ICEFREED.66
& rhosnow, ! IN SI density of snow ICEFREED.67
& cw ! IN quadratic drag coefficient ICEFREED.68
&,hicestop ! in max hice at which convergence is allowed. ICEFREED.69
&,hiceslow ! in min hice at which convergence is impeded. ICEFREED.70
REAL ICEFREED.71
& coriolis(imt,jmt) ! IN coriolis parameter ICEFREED.72
LOGICAL ICEFREED.73
& icy(imt,jmt), ! IN true if sea ice present ICEFREED.74
& ocean(imt,jmt) ! IN true if not land ICEFREED.75
REAL ICEFREED.76
& aice(imt,jmt), ! IN ice concentration ICEFREED.77
& hice(imt,jmt), ! IN GBM ice depth ICEFREED.78
& hsnow(imt,jmt), ! IN snow depth ICEFREED.79
& ucurrent(imt,jmtm1), ! IN x-cpt ocean surface current ICEFREED.80
& vcurrent(imt,jmtm1), ! IN y-cpt ocean surface current ICEFREED.81
& wsx_ice(imt,jmtm1), ! IN x-cpt GBM windstress on ice ICEFREED.82
& wsy_ice(imt,jmtm1), ! IN y-cpt GBM windstress on ice ICEFREED.83
& uice(imt,jmtm1), ! OUT x-cpt ice velocity ICEFREED.84
& vice(imt,jmtm1), ! OUT y-cpt ice velocity ICEFREED.85
& isx(imt,jmtm1), ! OUT x-cpt GBM ice stress on ocean ICEFREED.86
& isy(imt,jmtm1), ! OUT y-cpt GBM ice stress on ocean ICEFREED.87
& twaterx(imt,jmtm1), ! OUT x-cpt ocean to ice stress ICEFREED.88
& twatery(imt,jmtm1), ! OUT x-cpt ocean to ice stress ICEFREED.89
& mfx(imt,jmtm1), ! OUT x-cpt Coriolis stress ICEFREED.90
& mfy(imt,jmtm1), ! OUT y-cpt Coriolis stress ICEFREED.91
& inisx(imt,jmtm1), ! OUT x-cpt internal ice stress ICEFREED.92
& inisy(imt,jmtm1), ! OUT y-cpt internal ice stress ICEFREED.93
& radius, ! Radius of the earth ICEFREED.94
& dphi, ! Latitudinal grid space ICEFREED.95
& dlambda, ! Longitudinal grid space ICEFREED.96
& phit(jmt) ! Latitude on Tracer grid ICEFREED.97
! Declare local variables ICEFREED.98
INTEGER ICEFREED.99
& i,j ICEFREED.100
REAL ICEFREED.101
& rhow_cw, ICEFREED.102
& recip_rhow_cw, ICEFREED.103
& Tmag, ICEFREED.104
& r1, ICEFREED.105
& r2, ICEFREED.106
& cwstar, ICEFREED.107
& mf, ICEFREED.108
& c1, ! mass / unit area x coriolis param ICEFREED.109
& u1, ICEFREED.110
& v1, ICEFREED.111
& p ! Local variable to save ICEFREED.112
! repeated calc.s ICEFREED.113
LOGICAL ICEFREED.114
& icy_uv(imt,jmtm1) ICEFREED.115
REAL ICEFREED.116
& mass(imt,jmtm1), ! mass/area of ice on vely grid ICEFREED.117
& htrue(imt,jmt), ! actual ice depth ICEFREED.118
& htrue_uv(imt,jmtm1), ! actual ice depth on vely grid ICEFREED.119
& hsnow_uv(imt,jmtm1), ! snow depth on vely grid ICEFREED.120
& aice_uv(imt,jmtm1), ! ice conc. on vely grid ICEFREED.121
& wsx(imt,jmtm1), ! x-cpt windstress over ice ICEFREED.122
& wsy(imt,jmtm1), ! y-cpt windstress over ice ICEFREED.123
& Tx(imt,jmtm1), ICEFREED.124
& Ty(imt,jmtm1), ICEFREED.125
& u1mag(imt,jmtm1) ICEFREED.126
C ICEFREED.127
! 0 Set variables for calculations on velocity points ICEFREED.128
! 0.1 Initialise arrays ICEFREED.129
! 0.2 Calculate real ice depth ICEFREED.130
IF (J_JMT+J_OFFSET.EQ.JMT_GLOBAL) THEN ICEFREED.131
DO i=1,imt ICEFREED.132
htrue(i,j_jmt)=0.0 ICEFREED.133
ENDDO ICEFREED.134
ENDIF ICEFREED.135
DO J=J_1,J_JMT ICEFREED.136
DO i=1,imt ICEFREED.137
IF (icy(i,j)) then ICEFREED.138
htrue(i,j) = hice(i,j)/aice(i,j) ICEFREED.139
ELSE ICEFREED.140
htrue(i,j) = 0.0 ICEFREED.141
ENDIF ICEFREED.142
ENDDO ICEFREED.143
ENDDO ICEFREED.144
! 0.3 Interpolate prognostics to velocity points ICEFREED.145
CALL h_to_uv( ICEFREED.146
*CALL ARGOINDX ICEFREED.147
& htrue,htrue_uv,imt,jmt,jmtm1) ICEFREED.148
CALL h_to_uv( ICEFREED.149
*CALL ARGOINDX ICEFREED.150
& hsnow,hsnow_uv,imt,jmt,jmtm1) ICEFREED.151
CALL h_to_uv( ICEFREED.152
*CALL ARGOINDX ICEFREED.153
& aice ,aice_uv ,imt,jmt,jmtm1) ICEFREED.154
*IF DEF,MPP ICEFREED.155
! Ensure ICY and OCEAN contain values in halos. ICEFREED.156
CALL SWAPBOUNDS(ICY,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ICEFREED.157
CALL SWAPBOUNDS(OCEAN,IMT,JMT,O_EW_HALO,O_NS_HALO,1) ICEFREED.158
*ENDIF ICEFREED.159
! 0.4 'Set up' calculations on the velocity grid ICEFREED.160
DO J=J_1,J_JMTM1 ICEFREED.161
DO i=1,imtm1 ICEFREED.162
icy_uv(i,j) = ( ICEFREED.163
& icy(i,j).or.icy(i,j+1).or.icy(i+1,j).or.icy(i+1,j+1) ) ICEFREED.164
& .and.ocean(i,j).and.ocean(i+1,j).and.ocean(i,j+1) ICEFREED.165
& .and.ocean(i+1,j+1) ICEFREED.166
ENDDO ! i ICEFREED.167
IF (l_ocyclic) THEN ICEFREED.168
icy_uv(imt,j) = icy_uv(2,j) ICEFREED.169
ELSE ICEFREED.170
icy_uv(imt,j) = ( ICEFREED.171
& (icy(imt,j).or.icy(imt,j+1) ) ICEFREED.172
& .and.ocean(imt,j).and.ocean(imt,j+1) ) ICEFREED.173
ENDIF ICEFREED.174
ENDDO ! J ICEFREED.175
ICEFREED.176
rhow_cw=rhowater*cw ICEFREED.177
recip_rhow_cw=1.0/rhow_cw ICEFREED.178
DO J=J_1,J_JMTM1 ICEFREED.179
DO i=1,imt ICEFREED.180
IF (icy_uv(i,j)) THEN ICEFREED.181
mass(i,j) = rhoice*htrue_uv(i,j) ICEFREED.182
& + rhosnow*hsnow_uv(i,j) ICEFREED.183
wsx(i,j)=wsx_ice(i,j)/aice_uv(i,j) ICEFREED.184
wsy(i,j)=wsy_ice(i,j)/aice_uv(i,j) ICEFREED.185
! 0.5 RHS term (independent of ice velocity) ICEFREED.186
mf=mass(i,j)*coriolis(i,j) ICEFREED.187
Tx(i,j) = wsx(i,j)+mf*vcurrent(i,j) ICEFREED.188
Ty(i,j) = wsy(i,j)-mf*ucurrent(i,j) ICEFREED.189
! 1 Get Ice Velocity ICEFREED.190
! 1.1 Magnitude ICEFREED.191
p=mf*mf*recip_rhow_cw*recip_rhow_cw ICEFREED.192
u1mag(i,j) = SQRT( ICEFREED.193
& 0.5*(SQRT(p*p+ 4.0*recip_rhow_cw*recip_rhow_cw* ICEFREED.194
& (Tx(i,j)*Tx(i,j)+Ty(i,j)*Ty(i,j)))-p) ICEFREED.195
& ) ICEFREED.196
! 1.2 Components ICEFREED.197
cwstar = rhow_cw * u1mag(i,j) ICEFREED.198
c1 = 1.0/( cwstar*cwstar + mf*mf ) ICEFREED.199
u1 = c1 * ( cwstar*Tx(i,j) + mf *Ty(i,j) ) ICEFREED.200
v1 = c1 * ( -mf *Tx(i,j) + cwstar*Ty(i,j) ) ICEFREED.201
uice(i,j) = u1 + ucurrent(i,j) ICEFREED.202
vice(i,j) = v1 + vcurrent(i,j) ICEFREED.203
ELSE ICEFREED.204
uice(i,j)=0.0 ICEFREED.205
vice(i,j)=0.0 ICEFREED.206
ENDIF ICEFREED.207
ENDDO ICEFREED.208
ENDDO ICEFREED.209
call cnvstop( ICEFREED.210
*CALL ARGOINDX ICEFREED.211
& imt,imtm1,imtm2,jmt,jmtm1,jmtm2, ICEFREED.212
& uice,vice,hice,radius,dphi,dlambda ICEFREED.213
& ,phit,hicestop,hiceslow) ICEFREED.214
! 2 GBM Quadratic Ice-Ocean Stress ICEFREED.215
! Calculate stress exerted on the ocean by the ice ICEFREED.216
DO J=J_1,J_jmtm1 ICEFREED.217
DO i=1,imt ICEFREED.218
IF (icy_uv(i,j)) THEN ICEFREED.219
u1 = uice(i,j) - ucurrent(i,j) ICEFREED.220
v1 = vice(i,j) - vcurrent(i,j) ICEFREED.221
u1mag(i,j)=sqrt(u1*u1+v1*v1) ICEFREED.222
twaterx(i,j) = -rhow_cw * u1mag(i,j) * u1 ICEFREED.223
twatery(i,j) = -rhow_cw * u1mag(i,j) * v1 ICEFREED.224
isx(i,j) = -twaterx(i,j) * aice_uv(i,j) ICEFREED.225
isy(i,j) = -twatery(i,j) * aice_uv(i,j) ICEFREED.226
mf=mass(i,j)*coriolis(i,j) ICEFREED.227
mfx(i,j) = mf*vice(i,j) ICEFREED.228
mfy(i,j) = -mf*uice(i,j) ICEFREED.229
inisx(i,j) = -wsx(i,j)-twaterx(i,j)-mfx(i,j) ICEFREED.230
inisy(i,j) = -wsy(i,j)-twatery(i,j)-mfy(i,j) ICEFREED.231
ELSE ICEFREED.232
twaterx(i,j) = 0.0 ICEFREED.233
twatery(i,j) = 0.0 ICEFREED.234
isx(i,j) = 0.0 ICEFREED.235
isy(i,j) = 0.0 ICEFREED.236
mfx(i,j)=0.0 ICEFREED.237
mfy(i,j)=0.0 ICEFREED.238
inisx(i,j)=0.0 ICEFREED.239
inisy(i,j)=0.0 ICEFREED.240
ENDIF ICEFREED.241
ENDDO ICEFREED.242
ENDDO ICEFREED.243
*IF DEF,MPP ICEFREED.244
CALL SWAPBOUNDS(isx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.245
CALL SWAPBOUNDS(isy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.246
CALL SWAPBOUNDS(twaterx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.247
CALL SWAPBOUNDS(twatery,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.248
CALL SWAPBOUNDS(mfx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.249
CALL SWAPBOUNDS(mfy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.250
CALL SWAPBOUNDS(inisx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.251
CALL SWAPBOUNDS(inisy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1) ICEFREED.252
*ENDIF ICEFREED.253
RETURN ICEFREED.254
END ICEFREED.255
*ENDIF ICEFREED.256