*/
*/ This mod leans heavily upon Doug Cresswells EVP mod, but:
*/
*/ - it corrects wsx_ice and tmass to /aice[_uv]
*/ - it corrects isx output to *aice
*/ - it does the fix in freedrift-correct-uice_local.upd too, so don't include that!
*/
*/ Remember *not* to use "fixoi.upd" (or equivalent) to fix the problems with the
*/ ice-ocean stresses, since this version does that in this mod.
*/
*/ It also creates the variable "phi_t" which is to be used (if desired) to
*/ introduce the open-water creation param if Stern et al. See mod stern_c2d.upd
*/ This mod creates the variable within ice_ctl; you need stern_c2d to make use of it.
*/
*/ Creates ICE_EVP.f and TYPIEVP,ARGIEVP comdecks. Lower down we mod ICE_CTL.
*/ Even lower are various scattered mods to bring signe11 etc down as prognostics, etc.
*/
*/ You want to be using the userstashfile evp_vn3, not the combination of
*/ evp_vn2 and userpre-freedriftdiags2
*/
*/ Vn 1.1 - wmc 2003/02/28 - fix definition of phi_t (err, include the 22 terms :-)
*/
*COMDECK ARGIEVP
C
C ======================= COMDECK ARGIEVP =============================
C Arguments passed around within EVP module
c Ice model / UM wide variables
& L_OCYCLIC,L_ICESSTILT,L_ICYNPOL
&,IMT,IMTM1,IMTM2,JMT,JMTM1,J_1,J_2,J_JMT,J_JMTM1
&,O_EW_HALO,O_NS_HALO
&,ndte,dtts,dte,dter
&,t_dx,t_dxr,t_dy,t_dyr,cs,cst,csr,tng,radius_si,HT_E,HT_N
&,ocean,icy,icy_uv,icy_uvp,ocean_t_mask,ocean_uv_mask
&,uice,vice,aice,hice,hsnow,rhoice,rhosnow,rho_water_si
&,CORIOLIS,UCURRENT,VCURRENT
&,wsx,wsy,isx,isy
&,sstiltx,sstilty
&,waterx,watery
&,quad_ice_drag,
c Internal EVP variables
& ecc2,ecci,eccm,eccp,cw,sw,
& eyc, cstar, pstar,
& etan, etas, etae, etaw, prss, umass,fm,
& sig11ne, sig11se, sig11sw, sig11nw,
& sig12ne, sig12se, sig12sw, sig12nw,
& sig22ne, sig22se, sig22sw, sig22nw,
c Temporary EVP variables
& a2na, a2sa, a2ea, a2wa, b2n, b2s, b2e, b2w,
& t_dx8, t_dy8, edy, edx, eHN, eHE, eHNm, eHEm, HT_N4, HT_E4,
& h2n, h2s, h2e, h2w, prssn, prsss, prsse, prssw,
& Tdter,tdamp2,c2n,d2n,acoef,bcoef,ccoef,dcoef,ecoef,
& dxtr4,dytr4
C ===================== END OF COMDECK ARGIEVP =========================
*/
*COMDECK TYPIEVP
C
C ======================= COMDECK TYPIEVP =============================
LOGICAL
& L_OCYCLIC
&,L_ICESSTILT
&,L_ICYNPOL
INTEGER
& IMT,IMTM1,IMTM2 !No. points in row, minus 1, and minus 2
&,JMT,JMTM1 !No. Tracer (Velocity) rows on processor
&,J_1 !First non-halo row on processor
&,J_2
&,J_JMT !Last non-halo tracer row on processor
&,J_JMTM1 !Last non-halo velocity row on processor
&,O_EW_HALO,O_NS_HALO
REAL
& ndte !No. EVP subcycles per timestep
REAL
& dtts !Timestep of forcing
&,dte,dter !Subcycling timestep, and reciprocal
&,t_dx(IMT,JMT) !Grid lengths between tracer points, & recip.s
&,t_dxr(IMT,JMT) ,t_dy(IMT,JMT) ,t_dyr(IMT,JMT)
&,cs(JMTM1),cst(JMT),csr(jmtm1) !Cos(phi) on U and T grid
&,tng(jmt) !Tan(lat) U grid (dimension jmt?) see TYPOCONE
&,radius_si !Radius of the earth
&,HT_E(IMT,JMT) !Length of eastern edge of T-cell
&,HT_N(IMT,JMT) !Length of northern edge of T-cell
LOGICAL
& ocean(imt,jmt) !Land-sea mask (Tracer grid)
&,icy(imt,jmt) !Ice mask
&,icy_uv(imt,jmtm1)
&,icy_uvp(imt,jmtm1)
REAL
& ocean_t_mask(IMT,JMT)
&,ocean_uv_mask(IMT,JMTM1)
&,uice(IMT,JMTM1) ,vice(IMT,JMTM1)
&,aice(IMT,JMT) ,hice(IMT,JMT) ,hsnow(IMT,JMT)
&,rhoice,rhosnow,rho_water_si
&,coriolis(IMT,JMT)
&,ucurrent(IMT,JMTM1)
&,vcurrent(IMT,JMTM1)
&,wsx(IMT,JMTM1)
&,wsy(IMT,JMTM1)
&,isx(IMT,JMTM1)
&,isy(IMT,JMTM1)
&,sstiltx(IMT,JMTM1),sstilty(IMT,JMTM1)
&,waterx(IMT,JMTM1),watery(IMT,JMTM1)
&,quad_ice_drag
real ecc2 ! squared eccentricity of yield curve
real ecci ! 4./ecc2
real eccm ! 2.-0.5*ecci)
real eccp ! 1.+ecci*0.25
real cw ! cosine of turn. angle in water (phiwater)
real sw ! sine of turn. angle in water (phiwater)
real eyc ! coefficient for calc. the parameter E
real cstar ! c*, parameter in the defn of prss
real pstar ! P*, param in defn of prss (dyn/cm^2)
real etan(imt,jmt) ! shear viscosity in n triangle of T-cell
real etas(imt,jmt) ! shear viscosity in s triangle of T-cell
real etae(imt,jmt) ! shear viscosity in e triangle of T-cell
real etaw(imt,jmt) ! shear viscosity in w triangle of T-cell
real prss(imt,jmt) ! pressure P (centered in T-cell) (dyn/cm)
real umass(imt,jmtm1) ! mass of U-cell (kg)
real fm(imt,jmtm1) ! mass*coriolis param. on U grid
real sig11ne(imt,jmt) ! int. ice stress tensor, sigma_11, north
real sig11se(imt,jmt) ! sigma_11, east (g/s^2)
real sig11sw(imt,jmt) ! sigma_11, south
real sig11nw(imt,jmt) ! sigma_11, west
real sig12ne(imt,jmt) ! sigma_12, north
real sig12se(imt,jmt) ! sigma_12, east
real sig12sw(imt,jmt) ! sigma_12, south
real sig12nw(imt,jmt) ! sigma_12, west
real sig22ne(imt,jmt) ! sigma_22, north
real sig22se(imt,jmt) ! sigma_22, east
real sig22sw(imt,jmt) ! sigma_22, south
real sig22nw(imt,jmt) ! sigma_22 west
real
& a2na(imt,jmt), a2sa(imt,jmt), a2ea(imt,jmt), a2wa(imt,jmt),
& b2n(imt,jmt), b2s(imt,jmt), b2e(imt,jmt), b2w(imt,jmt),
& t_dx8(imt,jmt), t_dy8(imt,jmt),
& edy(imt,jmt), edx(imt,jmt), eHN(imt,jmt), eHE(imt,jmt),
& eHNm(imt,jmt), eHEm(imt,jmt), HT_N4(imt,jmt), HT_E4(imt,jmt),
& h2n(imt,jmt), h2s(imt,jmt), h2e(imt,jmt), h2w(imt,jmt),
& prssn(imt,jmt), prsss(imt,jmt), prsse(imt,jmt), prssw(imt,jmt)
real
& Tdter,tdamp2,c2n,d2n,acoef,bcoef,ccoef,dcoef,ecoef,
& dxtr4(imt,jmt),dytr4(imt,jmt)
integer i,j !local counters
C ===================== END OF COMDECK TYPIEVP ========================
*/
*/ Now the code...
*/
*DECK ICE_EVP
subroutine ICE_EVP(
*CALL ARGSTS
& L_OCYCLIC,L_ICESSTILT,L_ICYNPOL
&,IMT,JMT
&,J_1,J_2,J_JMT,J_JMTM1
&,O_EW_HALO,O_NS_HALO
&,dtts
&,dxt,dxtr,dyt,dytr,cs,cst,csr,tng,radius_si
&,ocean,icy
&,uice,vice
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
&,aice,hice,hsnow,rhoice,rhosnow,rho_water_si
&,CORIOLIS,UCURRENT,VCURRENT
&,wsx,wsy,isx,isy
&,sstiltx,sstilty
&,quad_ice_drag
&,dte,cw,sw,ecc2,eyc,cstar,pstar
&,STASHWORK
&,sw_len32
&,im_index
&,O_MYPE,O_NPROC
&,JFIN,JMT_GLOBAL
&,phi_t ! WMC stern (16) see stern_c2d.upd
& )
implicit none
*CALL CSUBMODL
*CALL PARPARM
*CALL TYPSIZE
*CALL TYPSTS
*CALL TYPIEVP
real phi_t(imt,jmt)
integer sw_len32,im_index
integer O_MYPE,O_NPROC
integer JFIN,JMT_GLOBAL
real STASHWORK(sw_len32)
real uorel,vorel,vrel
logical ocean_uv(imt,jmtm1)
real
& dsig11dx(IMT,JMTM1)
&,dsig12dy(IMT,JMTM1)
&,dsig21dx(IMT,JMTM1)
&,dsig22dy(IMT,JMTM1)
&,tmass(IMT,JMT)
C
&,uice_save(IMT,JMTM1)
&,vice_save(IMT,JMTM1)
C
real
& dxt(IMT),dxtr(IMT) !Grid lengths and reciprocals
&,dyt(JMT),dytr(JMT)
&,uice_init(IMT,JMTM1)
&,vice_init(IMT,JMTM1)
&,sstilt_stress_x(IMT,JMTM1)
&,sstilt_stress_y(IMT,JMTM1)
real coslamda(imt),sinlamda(imt)
integer k !sub-step loop counter
integer nlev !number of levels for stash
C Start of code
nlev=1
c CALL swap_1d(wsx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(wsy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(coriolis,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c Precalculate coslamda and sinlamda for polar island solution
do i=1,imt
coslamda(i)=cos(2*3.141592654*i/imtm2)
sinlamda(i)=sin(2*3.141592654*i/imtm2)
enddo
do j=1,JMTM1
do i=1,IMT
uice_init(i,j)=uice(i,j)
vice_init(i,j)=vice(i,j)
enddo
enddo
call evp_setup( ! set up var.s, hard-wire const.s init.ly
*CALL ARGIEVP
&,dxt,dxtr,dyt,dytr
&,O_MYPE,O_NPROC
&,OCEAN_UV
& )
call mass_prss( ! mass, pressure
*CALL ARGIEVP
&,tmass
& )
c Add and calc. diagnostic of sstilt stress, and sstilt
if (L_ICESSTILT) then
do j=1,jmtm1
do i=1,imt
sstilt_stress_x(i,j)=sstiltx(i,j)*umass(i,j)*9.81
sstilt_stress_y(i,j)=sstilty(i,j)*umass(i,j)*9.81
enddo
enddo
IF (SF(309,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,sstiltx,ocean,STASHWORK(SI(309,32,im_index)))
endif
IF (SF(310,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,sstilty,ocean,STASHWORK(SI(310,32,im_index)))
endif
IF (SF(311,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,sstilt_stress_x,icy,STASHWORK(SI(311,32,im_index)))
endif
IF (SF(312,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,sstilt_stress_y,icy,STASHWORK(SI(312,32,im_index)))
endif
endif
IF (SF(313,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,ucurrent,ocean,STASHWORK(SI(313,32,im_index)))
endif
IF (SF(314,32)) then
C WMC
C CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
CALL COPYODIAGL(IMT,JMTm1,nlev,.TRUE.,0.
& ,vcurrent,ocean,STASHWORK(SI(314,32,im_index)))
endif
IF (SF(300,32)) then
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
& ,tmass,icy,STASHWORK(SI(300,32,im_index)))
endif
IF (SF(301,32)) then
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
& ,prss,icy,STASHWORK(SI(301,32,im_index)))
endif
call stressprep( ! coefficients for stress calculations
*CALL ARGIEVP
& )
C ---------------------------------------------------------------
C EVP subcycling loop
C ---------------------------------------------------------------
DO k=1,ndte
call stress( ! stress tensor equation (sigma)
*CALL ARGIEVP
c WMC: not needed and not consistent!
c &,coslamda,sinlamda
c &,JFIN,JMT_GLOBAL
& ,phi_t
& )
call stepu( ! total surface stress, momentum eqn
*CALL ARGIEVP
&,uice_init,vice_init
&,dsig11dx
&,dsig12dy
&,dsig21dx
&,dsig22dy
&,JFIN,JMT_GLOBAL
&,coslamda,sinlamda
& )
C We may want the change-in-last-iteration diagnostics.
C Do this, of course, on second-to-last-loop only
IF (k .eq. ndte-1) THEN
IF (sf(326,32) .or. sf(328,32)) THEN
do j=2,J_JMTM1
do i=1,IMTM1
uice_save(i,j)=uice(i,j)
enddo
enddo
ENDIF
IF (sf(327,32) .or. sf(328,32)) THEN
do j=2,J_JMTM1
do i=1,IMTM1
vice_save(i,j)=vice(i,j)
enddo
enddo
ENDIF
ENDIF
END DO
C ---------------------------------------------------------------
C End of EVP subcycling loop
C ---------------------------------------------------------------
C Change-in-last-iteration diagnostics get saved here
IF (sf(326,32) .or. sf(328,32)) THEN
do j=2,J_JMTM1
do i=1,IMTM1
C Note: what we save as stash item is the difference of the current
C and last iteration: this should be small if all is well
uice_save(i,j)=uice_save(i,j)-uice(i,j)
enddo
enddo
ENDIF
IF (sf(327,32) .or. sf(328,32)) THEN
do j=2,J_JMTM1
do i=1,IMTM1
vice_save(i,j)=vice_save(i,j)-vice(i,j)
enddo
enddo
ENDIF
IF (sf(326,32)) THEN
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,uice_save,ocean_uv,STASHWORK(SI(326,32,im_index)))
ENDIF
IF (sf(327,32)) THEN
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,vice_save,ocean_uv,STASHWORK(SI(327,32,im_index)))
ENDIF
IF (sf(328,32)) THEN
do j=2,J_JMTM1
do i=1,IMTM1
C Note: in order to save one array, we overwrite uice_save with
C sqrt(udiff^2+vdiff^2)
uice_save(i,j)=sqrt(uice_save(i,j)**2+vice_save(i,j)**2)
enddo
enddo
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,uice_save,ocean_uv,STASHWORK(SI(328,32,im_index)))
ENDIF
c..Write full ice-ocean drag into isx and isy
do j=2,J_JMTM1
do i=1,IMTM1
if (icy_uv(i,j)) then
uorel = UCURRENT(i,j) - uice(i,j)
vorel = VCURRENT(i,j) - vice(i,j)
vrel = quad_ice_drag*rho_water_si*
& sqrt(uorel**2 + vorel**2)
C
C WMC. Note these are absolute. They will be *aice by ice_ctl
C on leaving EVP.
C
C WMC. Note that ICEFREED thinks that the ice-ocean stress has the
C sign indicated by:
C
C u1 = uice(i,j) - ucurrent(i,j)
C twaterx(i,j) = -rhow_cw * u1mag(i,j) * u1
C isx(i,j) = -twaterx(i,j) * aice_uv(i,j)
C ie
C isx = c * (u_ice - u_water)
C
C ICEFREED is correct. So lets fix EVP here. This means we
C can (must) drop mod fixoi.upd and we MUST NOT negate the
C stresses in calcdiff etc.
C
C isx(i,j) = vrel*(uorel*cw - vorel*sw)
C isy(i,j) = vrel*(uorel*sw + vorel*cw)
isx(i,j) = -vrel*(uorel*cw - vorel*sw)
isy(i,j) = -vrel*(uorel*sw + vorel*cw)
else
isx(i,j) = 0.
isy(i,j) = 0.
endif
enddo
isx(IMT,j)=isx(2,j)
isy(IMT,j)=isy(2,j)
enddo
c CALL swap_1d(isx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(isy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(dsig11dx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(dsig12dy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(dsig21dx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(dsig22dy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
IF (SF(304,32)) then
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,dsig11dx,ocean_uv,STASHWORK(SI(304,32,im_index)))
endif
IF (SF(305,32)) then
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,dsig12dy,ocean_uv,STASHWORK(SI(305,32,im_index)))
endif
IF (SF(306,32)) then
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,dsig21dx,ocean_uv,STASHWORK(SI(306,32,im_index)))
endif
IF (SF(307,32)) then
CALL COPYODIAGL(IMT,JMTM1,nlev,.TRUE.,0.
& ,dsig22dy,ocean_uv,STASHWORK(SI(307,32,im_index)))
endif
return
end
subroutine evp_setup( ! set up variables, hard-wire const.s initly
*CALL ARGIEVP
&,dxt,dxtr,dyt,dytr
&,O_MYPE,O_NPROC
&,OCEAN_UV
& )
implicit none
*CALL TYPIEVP
LOGICAL OCEAN_UV(imt,jmtm1) !Logical land sea mask-u pts
REAL
& dxt(IMT),
& dxtr(IMT),
& dyt(JMT),
& dytr(JMT)
&,g_over_f(jmtm1)
real
& s0 !dummy
integer O_MYPE,O_NPROC
*CALL C_MDI
C
logical first_time_called /.true./
if(O_MYPE .lt. (O_NPROC/2.)) then
sw=-sw
endif
dter=1.0/dte
ndte=dtts*dter
tdamp2=2.*eyc*dtts
C
IF (first_time_called) THEN
first_time_called=.false.
write(6,*)'WMC parameters from ICE_CTL'
write(6,*)'dte=',dte
write(6,*)'cw=',cw
write(6,*)'sw=',sw
write(6,*)'ecc2=',ecc2
write(6,*)'eyc=',eyc
write(6,*)'cstar=',cstar
write(6,*)'pstar=',pstar
write(6,*)'tdamp=',0.5*tdamp2
write(6,*)'ndte=',ndte
ENDIF
Tdter=tdamp2*dter
c2n=Tdter+ecc2
ecoef=0.5*(ecc2-1.)
d2n=c2n-ecoef
ccoef=1./c2n
dcoef=ccoef/(Tdter+1.)
bcoef=dcoef*d2n
acoef=dcoef*ecoef
if (ndte .ne. int(ndte)) then
write(6,*)'Warning EVP subcycle length doesnt cleanly
& divide ocean timestep'
write(6,*)'EVP (sub)timestep=',dte
write(6,*)'Ocean timestep=',dtts
endif
c Eccentricity variables
ecci=4./ecc2
eccm=2.-0.5*ecci
eccp=1.+ecci*0.25
c Grid variables (convert from CGS to SI)
do j=1,JMT
do i=1,IMT
t_dx(i,j)=dxt(i)*cst(j)*0.01
t_dxr(i,j)=1.0/t_dx(i,j)
dxtr4(i,j)=0.25*t_dxr(i,j)
t_dy(i,j)=dyt(j)*0.01
t_dyr(i,j)=1.0/t_dy(i,j)
dytr4(i,j)=0.25*t_dyr(i,j)
HT_E(i,j)=t_dy(i,j)
HT_N(i,j)=dxt(i)*cs(j)*0.01
enddo
enddo
C Land/Sea, Ice/Noice masking
do i=1,IMT
do j=1,JMT
ocean_t_mask(i,j)=0.
if (ocean(i,j)) then
ocean_t_mask(i,j)=1.
endif
enddo
enddo
do i=1,IMTM1
do j=1,J_JMTM1
ocean_uv_mask(i,j)=0.
if (ocean(i,j).and.ocean(i+1,j)
& .and.ocean(i,j+1).and.ocean(i+1,j+1)) then
ocean_uv_mask(i,j)=1.
else
endif
ocean_uv(i,j)=ocean(i,j).and.ocean(i+1,j).and.ocean(i,j+1)
& .and.ocean(i+1,j+1)
icy_uv(i,j) = (
& icy(i,j).or.icy(i,j+1).or.icy(i+1,j).or.icy(i+1,j+1) )
& .and.ocean_uv(i,j)
enddo
enddo
if(L_OCYCLIC) then
do j=J_1,J_JMTM1
ocean_uv_mask(IMT,j)=ocean_uv_mask(2,j)
ocean_uv_mask(1,j)=ocean_uv_mask(IMTM1,j)
c ocean_t_mask(IMT,j)=ocean_t_mask(2,j)
icy_uv(IMT,j)=icy_uv(2,j)
icy_uv(1,j)=icy_uv(IMTM1,j)
enddo
else
do j=J_1,J_JMTM1
ocean_uv_mask(IMT,j)=0.
if (ocean(i,j).and.ocean(i,j+1)) then
ocean_uv_mask(i,j)=1.
endif
icy_uv(IMT,j)= (
& icy(i,j).or.icy(i,j+1))
& .and.ocean(i,j).and.ocean(i,j+1)
enddo
endif
do j=2,jmtm1
do i=2,imt ! ice extent mask
! + points north and east of ice pack (U-cells)
icy_uvp(i,j) = (icy_uv(i-1,j) .or. icy_uv(i,j-1)
& .or. icy_uv(i-1,j-1) .or. icy_uv(i,j))
enddo
enddo
if(L_OCYCLIC) then
do j=J_1,J_JMTM1
icy_uvp(1,j)=icy_uvp(IMTM1,j)
enddo
else
do j=J_1,J_JMTM1
icy_uvp(1,j)= (icy_uv(i,j-1).or. icy_uv(i,j))
enddo
endif
c CALL swap_1d(icy_uv,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(icy_uvP,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(ocean_uv_mask,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(ocean_t_mask,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sstiltx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sstilty,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(UCURRENT,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(VCURRENT,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c Initialise stress tensor
if(sig11ne(1,J_1) .eq. RMDI) then
C write(6,*)'Wobble RMDI=',RMDI
do j=J_1,J_JMT
do i=1,imt
s0 = -0.5*pstar*hice(i,j)
& *exp(-cstar*(1.-aice(i,j)))
sig11ne(i,j) = s0
sig11se(i,j) = s0
sig11sw(i,j) = s0
sig11nw(i,j) = s0
sig22ne(i,j) = s0
sig22se(i,j) = s0
sig22sw(i,j) = s0
sig22nw(i,j) = s0
sig12ne(i,j) = 0.
sig12se(i,j) = 0.
sig12sw(i,j) = 0.
sig12nw(i,j) = 0.
enddo
enddo
do j=J_1,J_JMTM1
do i=1,IMT
uice(i,j) = ucurrent(i,j)
vice(i,j) = vcurrent(i,j)
enddo
enddo
endif
c CALL swap_1d(sig11ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(uice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(vice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(t_dx,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(t_dxr,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(t_dy,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(t_dyr,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(HT_E,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(HT_N,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
return
end
subroutine mass_prss(
*CALL ARGIEVP
&,tmass
& )
c.. Computes ice mass, pressure (strength)
implicit none
*CALL TYPIEVP
integer i, j
real tmass(imt,jmt)
c.. total mass of ice and snow, centered in T-cell
C
C WMC. Note correction to calculation of tmass
C
do j=1,jmt
do i=1,imt
if ((aice(i,j) .gt. 0.0) .and. (aice(i,j) .le. 1.0)) then
tmass(i,j) = rhoice*hice(i,j)/aice(i,j) ! k/m^2
& +rhosnow*hsnow(i,j)
else
tmass(i,j) = 0.0
endif
tmass(i,j) = tmass(i,j)*ocean_t_mask(i,j) ! mask
enddo
enddo
c CALL swap_1d(tmass,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c.. mass centered at velocity nodes (U-cells)
do j=1,J_JMTM1
do i=1,imtm1
umass(i,j) = 0.25*(tmass(i,j) + tmass(i+1,j)
& + tmass(i,j+1) + tmass(i+1,j+1)) !kg/m^2
enddo
enddo
if(L_OCYCLIC) then
do j=1,J_JMTM1
umass(IMT,j)=umass(2,j)
enddo
endif
c.. pressure P
do j=1,jmt
do i=1,imt
prss(i,j) = ocean_t_mask(i,j)*pstar*hice(i,j)
& *exp(-cstar*(1.-aice(i,j))) ! g/s^2
enddo
enddo
c CALL swap_1d(umass,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(prss,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
return
end
subroutine stressprep(
*CALL ARGIEVP
& )
c.. Computes quantities used in the subroutines stress and stepu.
c.. This subroutine is cryptic - I apologize - but it made the code
c.. faster by almost a factor of 2 on a Cray-YMP. Here we compute
c.. quantities needed in the stress tensor (sigma) and momentum (u)
c.. equations, but which do not change during the subcycling. Many of
c.. these variables are grouped in parties of four, one for each
c.. triangle (north, east, south, west).
implicit none
*CALL TYPIEVP
integer i, j
real Econst,ey,e2,en,ee,es,ew,zn,ze,zs,zw,c2n,c2e,c2s,c2w,
& c3n,c3e,c3s,c3w,d2n,d2e,d2s,d2w,a2n,a2e,a2s,a2w,prs,g
g=9.7976
do j=2,JMT
do i=1,imt
HT_N4(i,j) = 0.25/HT_N(i,j)
HT_E4(i,j) = 0.25/HT_E(i,j)
t_dx8(i,j) = 0.125*t_dxr(i,j)
t_dy8(i,j) = 0.125*t_dyr(i,j)
enddo
enddo
do j=J_1,J_JMTM1
do i=1,imt
fm(i,j) = CORIOLIS(i,j)*umass(i,j) ! Coriolis * mass
c.. for water stress
waterx(i,j) = ocean_uv_mask(i,j)*
& (UCURRENT(i,j)*cw - VCURRENT(i,j)*sw)
watery(i,j) = ocean_uv_mask(i,j)*
& (VCURRENT(i,j)*cw + UCURRENT(i,j)*sw)
enddo
enddo
C WMC. Note that modifying wsxy is NAUGHTY
if (L_ICESSTILT) then
do j=J_1,J_JMTM1
do i=1,imt
wsx(i,j) = ocean_uv_mask(i,j)*
& (wsx(i,j) - g*umass(i,j)*sstiltx(i,j))
wsy(i,j) = ocean_uv_mask(i,j)*
& (wsy(i,j) - g*umass(i,j)*sstilty(i,j))
enddo
enddo
else
do j=J_1,J_JMTM1
do i=1,imt
wsx(i,j) = ocean_uv_mask(i,j)*wsx(i,j)
wsy(i,j) = ocean_uv_mask(i,j)*wsy(i,j)
enddo
enddo
endif
c Swapbounds fm,waterx and watery for completeness.
c This shouldn't be strictly nec. but avoids passing garbage around.
c CALL swap_1d(fm,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(waterx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(watery,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
return
end
subroutine stress(
*CALL ARGIEVP
& ,phi_t
& )
c.. Calculates the internal stress components, sigma_ij, in the four
c.. triangles of each cell.
implicit none
*CALL TYPIEVP
C
real phi_t(imt,jmt)
C
integer i, j
real dun,dus,due,duw,dvn,dvs,dve,dvw,cc,dd,
& xi11ne,xi12ne,xi22ne,xi11se,xi12se,xi22se,
& xi11sw,xi12sw,xi22sw,xi11nw,xi12nw,xi22nw,
& c4ne,c4se,c4sw,c4nw,c5ne,c5se,c5sw,c5nw,
& deltane,deltase,deltasw,deltanw,
& prsdne,prsdse,prsdsw,prsdnw,
& prs2ne,prs2se,prs2sw,prs2nw
real rcon,radiusr,metric(imt,jmt)
real prs_sig(imt,jmt)
real xi11pol,xi12pol,xi22pol,sig11pol,sig12pol,sig22pol
real deltapol,prsspol,prsdpol,prs2pol,c4pol,c5pol
real sig11pol,sig12pol,sig22pol
real coslamda(imt),sinlamda(imt)
integer jfin,jmt_global
c CALL swap_1d(sig11ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
rcon=4e-9
radiusr=1.0/radius_si
do j=1,jmtm1
do i=2,imt
metric(i,j)=radiusr*0.25*(tng(j)*(vice(i-1,j)+vice(i,j)) +
& tng(j-1)*(vice(i-1,j-1)+vice(i,j-1)))
enddo
enddo
do j=2,jmt
do i=2,imt
if (ocean(i,j)) then
if (icy(i,j)) then
c Include metric terms for spherical polar coord.s
dun = (uice(i,j) - uice(i-1,j))*t_dxr(i,j)
dus = (uice(i,j-1) - uice(i-1,j-1))*t_dxr(i,j)
due = cst(j)*(uice(i,j)*csr(j) - uice(i,j-1)*csr(j-1))
& *t_dyr(i,j)
duw = cst(j)*(uice(i-1,j)*csr(j) - uice(i-1,j-1)*csr(j-1))
& *t_dyr(i,j)
dvn = (vice(i,j) - vice(i-1,j))*t_dxr(i,j)
dvs = (vice(i,j-1) - vice(i-1,j-1))*t_dxr(i,j)
dve = (vice(i,j) - vice(i,j-1))*t_dyr(i,j)
dvw = (vice(i-1,j) - vice(i-1,j-1))*t_dyr(i,j)
c Calc. strain rates on each (ne,se,sw,nw) triangle
xi11ne = dun - metric(i,j)
xi12ne = 0.5*(due + dvn)
xi22ne = dve
xi11se = dus - metric(i,j)
xi12se = 0.5*(due + dvs)
xi22se = dve
xi11sw = dus - metric(i,j)
xi12sw = 0.5*(duw + dvs)
xi22sw = dvw
xi11nw = dun - metric(i,j)
xi12nw = 0.5*(duw + dvn)
xi22nw = dvw
c.. Delta (in the denominator of zeta, eta) ! 1/s
deltane = sqrt( (xi11ne**2+xi22ne**2)*eccp
& + xi12ne**2*ecci + xi11ne*xi22ne*eccm)
deltasw = sqrt( (xi11sw**2+xi22sw**2)*eccp
& + xi12sw**2*ecci + xi11sw*xi22sw*eccm)
deltanw = sqrt( (xi11nw**2+xi22nw**2)*eccp
& + xi12nw**2*ecci + xi11nw*xi22nw*eccm)
deltase = sqrt( (xi11se**2+xi22se**2)*eccp
& + xi12se**2*ecci + xi11se*xi22se*eccm)
C
C See stern_c2d.upd for details. Form phi_t ish from eq (16)
phi_t(i,j)= 1/2. *
& (
& (1/4.*(deltane+deltanw+deltasw+deltase))
& -
& (1/4.*(xi11ne+xi11sw+xi11nw+xi11se
& +xi22ne+xi22sw+xi22nw+xi22se))
& )
C
c.. replacement pressure P/Delta
prsdne = prss(i,j)/max(rcon,deltane)
prsdsw = prss(i,j)/max(rcon,deltasw)
prsdnw = prss(i,j)/max(rcon,deltanw)
prsdse = prss(i,j)/max(rcon,deltase)
c.. P/2
prs2ne = 0.5*prsdne*deltane
prs2sw = 0.5*prsdsw*deltasw
prs2nw = 0.5*prsdnw*deltanw
prs2se = 0.5*prsdse*deltase
c.. save replacement pressure for principal stress calculation
prs_sig(i,j) = 2.*prs2ne ! match triangle with that used in s
if (deltane.le.rcon) prs_sig(i,j) = 0. ! masks zetamax stres
c.. cryptic stuff
c4ne = Tdter*sig11ne(i,j) + prsdne*xi11ne - prs2ne
c4sw = Tdter*sig11sw(i,j) + prsdsw*xi11sw - prs2sw
c4nw = Tdter*sig11nw(i,j) + prsdnw*xi11nw - prs2nw
c4se = Tdter*sig11se(i,j) + prsdse*xi11se - prs2se
c5ne = Tdter*sig22ne(i,j) + prsdne*xi22ne - prs2ne
c5sw = Tdter*sig22sw(i,j) + prsdsw*xi22sw - prs2sw
c5nw = Tdter*sig22nw(i,j) + prsdnw*xi22nw - prs2nw
c5se = Tdter*sig22se(i,j) + prsdse*xi22se - prs2se
c.. finally
sig11ne(i,j) = (acoef*c5ne + c4ne*bcoef)*
& ocean_t_mask(i,j)
sig11sw(i,j) = (acoef*c5sw + c4sw*bcoef)*
& ocean_t_mask(i,j)
sig11nw(i,j) = (acoef*c5nw + c4nw*bcoef)*
& ocean_t_mask(i,j)
sig11se(i,j) = (acoef*c5se + c4se*bcoef)*
& ocean_t_mask(i,j)
sig22ne(i,j) = (acoef*c4ne + c5ne*bcoef)*
& ocean_t_mask(i,j)
sig22sw(i,j) = (acoef*c4sw + c5sw*bcoef)*
& ocean_t_mask(i,j)
sig22nw(i,j) = (acoef*c4nw + c5nw*bcoef)*
& ocean_t_mask(i,j)
sig22se(i,j) = (acoef*c4se + c5se*bcoef)*
& ocean_t_mask(i,j)
sig12ne(i,j) = ((Tdter*sig12ne(i,j) + prsdne*xi12ne)*ccoef)*
& ocean_t_mask(i,j)
sig12sw(i,j) = ((Tdter*sig12sw(i,j) + prsdsw*xi12sw)*ccoef)*
& ocean_t_mask(i,j)
sig12nw(i,j) = ((Tdter*sig12nw(i,j) + prsdnw*xi12nw)*ccoef)*
& ocean_t_mask(i,j)
sig12se(i,j) = ((Tdter*sig12se(i,j) + prsdse*xi12se)*ccoef)*
& ocean_t_mask(i,j)
endif
endif
enddo
enddo
if(L_OCYCLIC) then
do j=2,jmt
sig11ne(1,j)=sig11ne(IMTM1,j)
sig11se(1,j)=sig11se(IMTM1,j)
sig11sw(1,j)=sig11sw(IMTM1,j)
sig11nw(1,j)=sig11nw(IMTM1,j)
sig22ne(1,j)=sig22ne(IMTM1,j)
sig22se(1,j)=sig22se(IMTM1,j)
sig22sw(1,j)=sig22sw(IMTM1,j)
sig22nw(1,j)=sig22nw(IMTM1,j)
sig12ne(1,j)=sig12ne(IMTM1,j)
sig12se(1,j)=sig12se(IMTM1,j)
sig12sw(1,j)=sig12sw(IMTM1,j)
sig12nw(1,j)=sig12nw(IMTM1,j)
enddo
else
do j=2,jmt
sig11ne(1,j)=0.
sig11se(1,j)=0.
sig11sw(1,j)=0.
sig11nw(1,j)=0.
sig22ne(1,j)=0.
sig22se(1,j)=0.
sig22sw(1,j)=0.
sig22nw(1,j)=0.
sig12ne(1,j)=0.
sig12se(1,j)=0.
sig12sw(1,j)=0.
sig12nw(1,j)=0.
enddo
endif
c CALL swap_1d(sig11ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig11nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig12nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22ne,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22se,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22sw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(sig22nw,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
return
end
subroutine stepu(
*CALL ARGIEVP
&,uice_init,vice_init
&,dsig11dx
&,dsig12dy
&,dsig21dx
&,dsig22dy
&,JFIN,JMT_GLOBAL
&,coslamda,sinlamda
& )
c.. Calculation of the surface stresses
c.. Integration of the momentum equation to find velocity (u,v)
implicit none
*CALL TYPIEVP
integer i, j,iloc
real uorel,vorel,vrel,taux,tauy,umassdter,cca,ccb,ab2,
& s11,s12,s21,s22,c1,c2,vrel_v(imt*jmt)
real s11ns(imt,jmt-1), s11ew(imt,jmt)
real s22ew(imt-1,jmt), s22ns(imt,jmt)
real s12ewi(imt-1,jmt), s12ns(imt,jmt)
real s12nsj(imt,jmt-1), s12ew(imt,jmt)
real
& uice_init(IMT,JMTM1)
&,vice_init(IMT,JMTM1)
&,uice_pol_loc_cart !Locally cartesian ice velocity at north pole
&,vice_pol_loc_cart !
&,uice_pol_sph !Vorticity correction for polar row
&,vice_pol_sph !Divergence correction for polar row
real coslamda(imt),sinlamda(imt) !cos and sin of longitude
real
& dsig11dx(IMT,JMTM1)
&,dsig12dy(IMT,JMTM1)
&,dsig21dx(IMT,JMTM1)
&,dsig22dy(IMT,JMTM1)
integer JFIN,JMT_GLOBAL
real
& s1122(8,imt,jmt)
c!!!!!!!!!!!!COPY WSX INTO POLAR ROW _ TRIAL ONLY: NOT A FIX
if ( JFIN .ge. JMT_GLOBAL) then
do i=1,imt
wsx(i,j_jmt)=wsx(i,j_jmt-1)
wsy(i,j_jmt)=wsy(i,j_jmt-1)
enddo
endif
c.. "some cryptic but useful arrays"-E.Hunke
do j=1,jmt
do i=1,imt
if (icy(i,j)) then
s1122(1,i,j) = (sig11nw(i,j) + sig11ne(i,j))*dxtr4(i,j)
s1122(2,i,j) = (sig11sw(i,j) + sig11se(i,j))*dxtr4(i,j)
s1122(3,i,j) = (sig22se(i,j) + sig22ne(i,j))*dytr4(i,j)
s1122(4,i,j) = (sig22sw(i,j) + sig22nw(i,j))*dytr4(i,j)
s1122(5,i,j) = (sig12nw(i,j) + sig12ne(i,j))*dxtr4(i,j)
s1122(6,i,j) = (sig12sw(i,j) + sig12se(i,j))*dxtr4(i,j)
s1122(7,i,j) = (sig12se(i,j) + sig12ne(i,j))*dytr4(i,j)
s1122(8,i,j) = (sig12sw(i,j) + sig12nw(i,j))*dytr4(i,j)
else
s1122(1,i,j) = 0.
s1122(2,i,j) = 0.
s1122(3,i,j) = 0.
s1122(4,i,j) = 0.
s1122(5,i,j) = 0.
s1122(6,i,j) = 0.
s1122(7,i,j) = 0.
s1122(8,i,j) = 0.
endif
enddo
enddo
c.. integrate the momentum equation
iloc=0
do j=2,J_JMTM1
do i=1,IMTM1
if (icy_uv(i,j)) then
uorel = UCURRENT(i,j) - uice(i,j)
vorel = VCURRENT(i,j) - vice(i,j)
iloc=iloc+1
C WMC vrel_v(iloc) = uorel**2 + vorel**2 !
vrel_v(iloc) = sqrt(uorel**2 + vorel**2)
endif
enddo
enddo
C WMC call sqrt_v(iloc,vrel_v,vrel_v)
iloc=0
do j=2,J_JMTM1
do i=1,IMTM1
if(icy_uv(i,j)) then
iloc=iloc+1
vrel=quad_ice_drag*rho_water_si*vrel_v(iloc)
c.. ice/ocean stress
isx(i,j) = vrel*waterx(i,j)
isy(i,j) = vrel*watery(i,j)
c.. TOTAL surface stress--includes wind stress, ice/ocean stress, tilt
C wmc
taux = wsx(i,j) + isx(i,j)
tauy = wsy(i,j) + isy(i,j)
c.. alpha, beta are defined in Hunke and Dukowicz (1997), section 3.2
umassdter = umass(i,j)*dter ! m/dte, g/cm^2 s
cca = umassdter + vrel*cw ! alpha, g/cm^2 s
ccb = fm(i,j) + vrel*sw ! beta, g/cm^2 s
ab2 = max(cca**2 + ccb**2,1.e-11)
c.. "more cryptic stuff"-E.Hunke
s11=s1122(2,i+1,j+1)-s1122(2,i,j+1)+s1122(1,i+1,j)-s1122(1,i,j)
s22=s1122(4,i+1,j+1)-s1122(4,i+1,j)+s1122(3,i,j+1)-s1122(3,i,j)
s21=s1122(6,i+1,j+1)-s1122(6,i,j+1)+s1122(5,i+1,j)-s1122(5,i,j)
s12=s1122(8,i+1,j+1)-s1122(8,i+1,j)+s1122(7,i,j+1)-s1122(7,i,j)
c.. finally, the velocity components
c1 = s11 + s12 + taux + umassdter*uice(i,j)
c2 = s21 + s22 + tauy + umassdter*vice(i,j)
uice(i,j) = (cca*c1 + ccb*c2)*ocean_uv_mask(i,j)/ab2 !
vice(i,j) = (cca*c2 - ccb*c1)*ocean_uv_mask(i,j)/ab2
else
c.. set velocity and stress to zero on land and (nearly) open water
uice(i,j) = 0.
vice(i,j) = 0.
isx(i,j) = 0.
isy(i,j) = 0.
s11=0.
s22=0.
s12=0.
s21=0.
endif
dsig11dx(i,j)=s11
dsig12dy(i,j)=s12
dsig21dx(i,j)=s21
dsig22dy(i,j)=s22
enddo
if (L_OCYCLIC) then
uice(IMT,j)=uice(2,j)
vice(IMT,j)=vice(2,j)
isx(IMT,j)=isx(2,j)
isy(IMT,j)=isy(2,j)
dsig11dx(IMT,j)=dsig11dx(2,j)
dsig12dy(IMT,j)=dsig12dy(2,j)
dsig21dx(IMT,j)=dsig21dx(2,j)
dsig22dy(IMT,j)=dsig22dy(2,j)
endif
enddo
if((JFIN .ge. JMT_GLOBAL).and.L_ICYNPOL) then
uice_pol_loc_cart=0.
vice_pol_loc_cart=0.
uice_pol_sph=0.
vice_pol_sph=0.
do i=1,imtm2
uice_pol_loc_cart=uice_pol_loc_cart
& + coslamda(i)*uice(i,j_jmtm1-1)
& - sinlamda(i)*vice(i,j_jmtm1-1)
vice_pol_loc_cart=vice_pol_loc_cart
& + sinlamda(i)*uice(i,j_jmtm1-1)
& + coslamda(i)*vice(i,j_jmtm1-1)
uice_pol_sph=uice_pol_sph+uice(i,j_jmtm1-1)
vice_pol_sph=vice_pol_sph+vice(i,j_jmtm1-1)
enddo
uice_pol_loc_cart=uice_pol_loc_cart/imtm2
vice_pol_loc_cart=vice_pol_loc_cart/imtm2
uice_pol_sph=uice_pol_sph/(3*imtm2)
vice_pol_sph=vice_pol_sph/(3*imtm2)
do i=1,imt
uice(i,j_jmtm1)=coslamda(i)*uice_pol_loc_cart
& + sinlamda(i)*vice_pol_loc_cart
& + uice_pol_sph
vice(i,j_jmtm1)=- sinlamda(i)*uice_pol_loc_cart
& + coslamda(i)*vice_pol_loc_cart
& + vice_pol_sph
enddo
endif
c CALL swap_1d(uice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(vice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(isx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(isy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
return
end
*/---------------------------------------------------------------------
*IDENT EVP
*/
*DECLARE ICE_CTL
*/
*/ 1. Fix uice/uice_local confusion
*/
*D ODC1F405.269,ODC1F405.272
C! Local variables for ice dynamics
C
C REAL
C & uice_local(imt_ice,jmtm1_ice)
C &,vice_local(imt_ice,jmtm1_ice)
*/
*/ 2. That done, now add in the EVP
*/
*I JT161193.153
real uice_fd(IMT,JMTM1) ! To allow us to compare to FD
real vice_fd(IMT,JMTM1) ! ditto, V
real uice_evp_fd(IMT,JMTM1) ! EVP-FD, U
real vice_evp_fd(IMT,JMTM1) ! EVP-FD, V
real wsx (IMT,JMTM1) ! wsx_ice/aice
real wsy (IMT,JMTM1) ! ditto y
real aice_uv (IMT,JMTM1) ! aice on U grid
real phi_t (IMT,JMT) ! extra open water source term
*/ From Dougs modset
*I OJC3F400.4
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
*I JT161193.121
&,sstiltx,sstilty
*I ODC1F405.257
REAL
& sig11ne(IMT,JMT) !Ice internal stresses
&,sig11se(IMT,JMT) !for EVP ice dynamics
&,sig11sw(IMT,JMT)
&,sig11nw(IMT,JMT)
&,sig12ne(IMT,JMT)
&,sig12se(IMT,JMT)
&,sig12sw(IMT,JMT)
&,sig12nw(IMT,JMT)
&,sig22ne(IMT,JMT)
&,sig22se(IMT,JMT)
&,sig22sw(IMT,JMT)
&,sig22nw(IMT,JMT)
*I ODC1F405.267
&,sstiltx(IMT,JMTM1),sstilty(IMT,JMTM1)
*D ODC1F405.254,ODC1F405.255
C Corrects dimensioning of WSXY_ICE
&,WSX_ICE(imt_idr,jmtm1_idr) ! IN ZONAL COMPONENT OF WIND STRESS.
&,WSY_ICE(imt_idr,jmtm1_idr) ! IN MERIDIONAL COMPONENT OF " "
*D ODC1F405.287,ODC1F405.340
*/
*/ Ice_ctl is such a mess: lets re-write the dynamics bit
*/
C We need to construct aice_uv. Code needs a tidy: icefreed does the
C same thing: really it should be done in ice_ctl only and once.
CALL h_to_uv(
*CALL ARGOINDX
& aice ,aice_uv ,imt,jmt,jmtm1)
DO J = J_1, J_JMTm1
DO I = 1, IMT
C Ughh.
if (aice_uv(i,j) .gt. 0.0 .and. aice_uv(i,j) .le. 1.0) then
wsx(i,j)=wsx_ice(i,j)/aice_uv(i,j)
wsy(i,j)=wsy_ice(i,j)/aice_uv(i,j)
else
wsx(i,j)=0.0
wsy(i,j)=0.0
endif
ENDDO
ENDDO
IF (L_ICEFREEDR .or. L_ICEEVP) THEN
IF (L_OTIMER) CALL TIMER('ICEFREEDR',3)
C Skip call to icefreed. For unknown reasons, it often crashes
C (fpe 8) under fujitsu, especially (but not only) when combined with ridging.
C Since this is only a diag, its not worth the pain.
GOTO 1234
CALL ICEFREEDR(
*CALL ARGOINDX
& imt,imtm1,imtm2,jmt,jmtm1,jmtm2,
& rhoice,rho_water_si,rhosnow,
& quad_ice_drag,hicestop,hiceslow,
& coriolis,
& icy,ocean,aice,hice,hsnow,
& ucurrent,vcurrent,
& wsx_ice,wsy_ice,
& uice_fd ,vice_fd ,
& isx,isy,
& twaterx,twatery,
& mfx,mfy,
& inisx,inisy,
& radius,dphi,dlambda,phit
& )
C
C We may want to save FD's isxy for comparison
C with EVP's
C
IF (SF(322,32)) CALL COPYODIAGL(IMT,JMTm1,1,.TRUE.,0.
& ,isx,ocean,STASHWORK(SI(322,32,im_index)))
IF (SF(323,32)) CALL COPYODIAGL(IMT,JMTm1,1,.TRUE.,0.
& ,isy,ocean,STASHWORK(SI(323,32,im_index)))
C
1234 CONTINUE
IF (L_OTIMER) CALL TIMER('ICEFREEDR',4)
ENDIF ! L_ICEFREEDR = true
IF (L_OTIMER) CALL TIMER('ICE_EVP',3)
CALL ICE_EVP(
*CALL ARGSTS
& L_OCYCLIC,L_ICESSTILT,L_ICYNPOL
&,IMT,JMT
&,J_1,J_2,J_JMT,J_JMTM1
&,O_EW_HALO,O_NS_HALO
&,dtts
&,dxt,dxtr,dyt,dytr,cs,cst,csr,tng,radius_si
&,ocean,icy
C WMC. Note: uvice, not uvice_evp
&,uice,vice
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
&,aice,hice,hsnow,rhoice,rhosnow,rho_water_si
&,CORIOLIS,UCURRENT,VCURRENT
c &,wsx_ice,wsy_ice,isx,isy
c WMC. Note call with wsx not wsx_ice
&,wsx,wsy,isx,isy
&,sstiltx,sstilty
&,quad_ice_drag
&,dte,cw,sw,ecc2,eyc,cstar,pstar
&,STASHWORK
&,sw_len32
&,im_index
&,O_MYPE,O_NPROC
&,JFIN,JMT_GLOBAL
&,phi_t ! WMC stern (16) see stern_c2d.upd
&)
IF (L_OTIMER) CALL TIMER('ICE_EVP',4)
C
C WMC. Convert EVP's isxy (absolute) to over-seaice-fraction
C
DO J = J_1, J_JMTm1
DO I = 1, IMT
isx(i,j)=isx(i,j)*aice_uv(i,j)
isy(i,j)=isy(i,j)*aice_uv(i,j)
ENDDO
ENDDO
C