head 1.2;
access;
symbols;
locks
wmc:1.2; strict;
comment @# @;
1.2
date 2002.07.01.13.11.45; author wmc; state Exp;
branches;
next 1.1;
1.1
date 2002.07.01.13.10.26; author wmc; state Exp;
branches;
next ;
desc
@Vaguely, the original
..
@
1.2
log
@Fix jmtm1 bugs with COPYODIAG
@
text
@*IDENT bench_EVP
*/
*/ Try to get ICE_SIMPLE/EVP to go together in the same piece of
*/ code. Main changes are some extra checks and an extra call
*/ involving icedrift...
*/
*/ WMC: minor changes.
*/ 2002/04/29: put sqrt_v inline (since sqrt_v doesn't exist...)
*/ 2002/04/29: comment out calls to swap_1d: we can';t use the cray version.
*/ on 2 procs or less, this is safe. 4 probably ok too.
*/ 2002/04/29: add diags for comparison against freedrift
*/ 2002/07/01: fix bug with jmt/jmtm1 in copyodiag that crashed 32-bit code
*/---------------------------------------------------------------------
*DECLARE BLOKCALC
*D ORH3F403.245
&,SURFTEMP,SURFSAL,NEWICE,UCURRENT,VCURRENT,UICE,VICE
*I ODC1F405.11
&,UICE(IMT_IDR,JMTM1_IDR)
&,VICE(IMT_IDR,JMTM1_IDR)
*D ORH5F401.6
&,AICE,AICE_UV,HICE,HSNOW,HICE_REF,CARYHEAT,ICY,FLXTOICE
*D ORH5F401.10
&,SURFTEMP, SURFSAL, NEWICE, UCURRENT, VCURRENT,UICE,VICE
*/---------------------------------------------------------------------
*DECLARE BLOKCNTL
*D ORH5F401.15
&,AICE,AICE_UV,HICE,HSNOW,HICE_REF,CARYHEAT,ICY,FLXTOICE
*D ORH5F401.19
&,SURFTEMP, SURFSAL, NEWICE, UCURRENT, VCURRENT,UICE,VICE
*I ORH1F305.3408
&, AICE_UV(*)
*I ORH1F305.3420
&,UICE(*)
&,VICE(*)
*D ORH5F401.25
&,AICE,AICE_UV,HICE,HSNOW,HICE_REF,CARYHEAT,ICY,FLXTOICE
*D ORH5F401.29
&,SURFTEMP, SURFSAL, NEWICE, UCURRENT, VCURRENT,UICE,VICE
*/---------------------------------------------------------------------
*DECLARE CLINIC
*D OOM1F405.749
&,FLXTOICE,FLXTOICEP
&,uice,vice,aice_uv
& )
*I OOM1F405.769
&,uice(imt_ice),vice(imt_ice),aice_uv(imt_ice)
*/
*/ modify calls to vdifcalc to accommodate both
*/ SIMPLE and EVP...
*/
*D ODC1F405.157
*D JT161193.364
& WSX,WSY,GM,gnu
&,uice,vice,aice_uv,quad_ice_drag,rho_water_si
&,L_ICEFREEDR,L_ICEEVP
&,wsx_leads,wsy_leads
& )
*D ORH1F305.3150,ORH1F305.3155
*/
*/---------------------------------------------------------------------
*DECLARE DERVSIZE
*I ORH1F305.4655
IF (L_ICEEVP) THEN
IMT_EVP=IMT
JMT_EVP=JMT
JMTM1_EVP=JMTM1
ELSE
IMT_EVP=1
JMT_EVP=1
JMTM1_EVP=1
ENDIF
*/---------------------------------------------------------------------
*DECLARE INITOI1A
*/
*D ODC1F405.344,ODC1F405.345
REAL wsx(imt,jmtm1)
REAL wsy(imt,jmtm1)
*D INITOI1A.11
&,aice,hice,newice )
*D INITOI1A.67
REAL aice(imt,jmt),hice(imt,jmt) !Ice conc/thickness
*I INITOI1A.111
hice(i,j)=0.
*I INITOI1A.122
hice(i,j)=0.
*I ORH7F402.85
&,fkmq
*I INITOI1A.57
REAL fkmq(imt,jmtm1) !Number of levels on u-v grid(land=0.0)
*I INITOI1A.148
do j=1,jmtm1
do i=1,imt
if(fkmq(i,j) .le. 0.1) then
wsx(i,j)=0.0
wsy(i,j)=0.0
endif
enddo
enddo
*I INITOI1A.165
c CALL swap_1d(wsx_ice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(wsy_ice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(wsx_leads,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(wsy_leads,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
*/---------------------------------------------------------------------
*DECLARE OCNFRST1
*I JT101193.9
&,D1(joc_isig11ne)
&,D1(joc_isig11se)
&,D1(joc_isig11sw)
&,D1(joc_isig11nw)
&,D1(joc_isig12ne)
&,D1(joc_isig12se)
&,D1(joc_isig12sw)
&,D1(joc_isig12nw)
&,D1(joc_isig22ne)
&,D1(joc_isig22se)
&,D1(joc_isig22sw)
&,D1(joc_isig22nw)
*/---------------------------------------------------------------------
*DECLARE OCN_CTL
*I ODC1F405.372
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
*/ *I ORH7F402.74
*B OCN_CTL.143
*CALL UMSCALAR
*D ODC1F405.375
&,WSX(IMT,JMTM1),WSY(IMT,JMTM1)
*I ODC1F405.374
REAL
& sig11ne(IMT_EVP,JMT_EVP) !Ice internal stresses
&,sig11se(IMT_EVP,JMT_EVP) !for EVP ice dynamics
&,sig11sw(IMT_EVP,JMT_EVP)
&,sig11nw(IMT_EVP,JMT_EVP)
&,sig12ne(IMT_EVP,JMT_EVP)
&,sig12se(IMT_EVP,JMT_EVP)
&,sig12sw(IMT_EVP,JMT_EVP)
&,sig12nw(IMT_EVP,JMT_EVP)
&,sig22ne(IMT_EVP,JMT_EVP)
&,sig22se(IMT_EVP,JMT_EVP)
&,sig22sw(IMT_EVP,JMT_EVP)
&,sig22nw(IMT_EVP,JMT_EVP)
*I ORH1F305.5551
&,sstiltx(IMT_ICE,JMTM1_ICE),sstilty(IMT_ICE,JMTM1_ICE)
&,d2pdtdy(IMT,JMTM1),d2pdtdx(IMT,JMTM1)
&,dtsfr,gravr,rhor
*I ORH7F402.77
&,fkmq
*I ORH1F305.5598
&,UICE,VICE
*I ORH1F305.5616
if(L_ICEEVP .and. L_ICESSTILT) then
dtsfr=1./dtsf
gravr=1./grav
rhor=1./RHOWATER
do j=1,J_JMTM1
do i=1,imtm1
d2pdtdy(i,j)=(0.5*(p(i,j+1)+p(i+1,j+1)-p(i,j)-p(i+1,j))
& -0.5*(pb(i,j+1)+pb(i+1,j+1)-pb(i,j)-pb(i+1,j)))
& *hr(i,j)*dyur(j)*dtsfr
d2pdtdx(i,j)=(0.5*(p(i+1,j)+p(i+1,j+1)-p(i,j)-p(i,j+1))
& -0.5*(pb(i+1,j)+pb(i+1,j+1)-pb(i,j)-pb(i,j+1)))
& *hr(i,j)*dxur(j)*dtsfr
sstiltx(i,j)=0.01*gravr*(0.5*dtsfr*zu(i,j)+d2pdtdy(i,j))
sstilty(i,j)=0.01*gravr*(0.5*dtsfr*zv(i,j)-d2pdtdx(i,j))
enddo
sstiltx(imt,j)=sstiltx(2,j)
sstilty(imt,j)=sstilty(2,j)
d2pdtdx(imt,j)=d2pdtdx(2,j)
d2pdtdy(imt,j)=d2pdtdy(2,j)
enddo
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(d2pdtdx,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(d2pdtdy,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
endif
*I JT161193.63
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
! sigmas - internal ice stresses for EVP ice dynamics
*I JT161193.72
&,sstiltx,sstilty
*D OJC2F400.81
&,aice,hice,newice )
*/---------------------------------------------------------------------
*DECLARE READNLST
*I ODC1F405.171
NAMELIST /ICEEVPNL/
& dte,cw,sw,ecc2,eyc,cstar,pstar
*I OCC0F400.53
IF (L_ICEEVP) then
READ(5,ICEEVPNL)
WRITE(6,ICEEVPNL)
ELSE
dte=0.0
cw=0.0
sw=0.0
ecc2=0.0
eyc=0.0
cstar=0.0
pstar=0.0
ENDIF
*/---------------------------------------------------------------------
*DECLARE ROWCALC
*D ORH5F401.32
&,AICE,AICE_UV,HICE,HSNOW,HICE_REF,CARYHEAT,ICY,FLXTOICE
*D ORH5F401.36
&,SURFTEMP, SURFSAL, NEWICE, UCURRENT, VCURRENT,UICE,VICE
*I ODC1F405.389
&,UICE(IMT_EVP,JMTM1_EVP)
&,VICE(IMT_EVP,JMTM1_EVP)
*D OOM1F405.690
&,FLXTOICE(1,J_ICE),FLXTOICE(1,J_ICEP1)
&,uice(1,J_ICE),vice(1,J_ICE),aice_uv(1,J_ICE)
&)
*/---------------------------------------------------------------------
*DECLARE ROW_CTL
*D OOM1F405.383
&,SURFTEMP,SURFSAL,NEWICE,UCURRENT,VCURRENT,UICE,VICE
&,ZTD,XF,YF,SWZVRT
*I ODC1F405.420
&,UICE(IMT_IDR,JMTM1_IDR)
&,VICE(IMT_IDR,JMTM1_IDR)
*D ORH3F403.254
&,SURFTEMP,SURFSAL,NEWICE,UCURRENT,VCURRENT,UICE,VICE
*/---------------------------------------------------------------------
*DECLARE STOCNPT1
*I GRB4F305.487
joc_isig11ne=SI(210,0,im_index)
joc_isig11se=SI(211,0,im_index)
joc_isig11sw=SI(212,0,im_index)
joc_isig11nw=SI(213,0,im_index)
joc_isig12ne=SI(214,0,im_index)
joc_isig12se=SI(215,0,im_index)
joc_isig12sw=SI(216,0,im_index)
joc_isig12nw=SI(217,0,im_index)
joc_isig22ne=SI(218,0,im_index)
joc_isig22se=SI(219,0,im_index)
joc_isig22sw=SI(220,0,im_index)
joc_isig22nw=SI(221,0,im_index)
*/---------------------------------------------------------------------
*DECLARE VDIFCALC
*I VDIFCALC.47
& ,uice,vice,aice_uv,quad_ice_drag,rho_water_si
& ,L_ICEFREEDR,L_ICEEVP
& ,wsx_leads,wsy_leads
*I VDIFCALC.62
REAL wsx_leads(imt),wsy_leads(imt)
REAL uice(imt),vice(imt),aice_uv(imt),quad_ice_drag,rho_water_si
LOGICAL L_ICEFREEDR,L_ICEEVP
*I VDIFCALC.78
&,cd_rho_ocn_ice_rel_sp(imt)
*I VDIFCALC.87
if(L_ICEFREEDR.or.L_ICEEVP) then
C Calculate relative velocity for inclusion of
C semi-implicit ocean-ice stress
C (also multipying by the drag coeffecient and density)
conv=10.
do i=1,imt
cd_rho_ocn_ice_rel_sp(i)=quad_ice_drag*rho_water_si
& *conv*aice_uv(i)
& *sqrt((uice(i)-ub(i,1))**2+(vice(i)-vb(i,1))**2)
enddo
endif
*I VDIFCALC.139
if(L_ICEFREEDR.or.L_ICEEVP) then
C Add ice-ocean stress terms in semi-implicitly
C (out of the quadratic term, the relative
C speed in explicit, and the velocity implicit)
do i=2,imtm1
efdr(i)=fxb/(bb(i,1)+(C2DTUV*cd_rho_ocn_ice_rel_sp(i)/DZ(1))
& -cc(i,1)*ee(i,2))
enddo
endif
*I VDIFCALC.154
if(L_ICEFREEDR.or.L_ICEEVP) then
C Add ice-ocean stress terms in semi-implicitly
C (out of the quadratic term, the relative
C speed in explicit, and the velocity implicit)
IF(M.EQ.1) THEN
DO 303 I=2,IMTM1
bb(I,1)=((WSX_leads(I) + cd_rho_ocn_ice_rel_sp(i)*uice(i))
& *conv*C2DTUV/DZ(1)+UB(I,1)+UA(I,1)*C2DTUV
& +cc(I,1)*ff(I,2,1))*efdr(I)
303 CONTINUE
ELSE
DO 306 I=2,IMTM1
bb(I,1)=((WSY_leads(I) + cd_rho_ocn_ice_rel_sp(i)*vice(i))
& *conv*C2DTUV/DZ(1)+VB(I,1)+VA(I,1)*C2DTUV
& +cc(I,1)*ff(I,2,2))*efdr(I)
306 CONTINUE
ENDIF
ELSE
*I VDIFCALC.165
ENDIF
*/---------------------------------------------------------------------
*/ADVECTION OVER THE POLAR ISLAND
*DECLARE ICEADVEC
*I ICEADVEC.47
&,L_ICYNPOL
*I ICEADVEC.73
&,L_ICYNPOL
*I ICEADVEC.117
if ((JFIN .ge. JMT_GLOBAL).and.L_ICYNPOL) then
c Remove polar island from land mask
do i=1,imt
hmask(i,j_jmt)=1.
hmask(i,jmt)=1.
enddo
endif
*I ICEADVEC.89
C Variables introduced for treatment of polar island
real imtm2r,total_jmt_hice,total_jmt_aice,total_jmt_hsnow
*// Remove other polar row code (retrun later with logical switch,...
*// ..., FOAM use it!)
*D ORH9F402.266,ICEADVEC.234
if (L_ICYNPOL) then
total_jmt_hice=0.
total_jmt_aice=0.
total_jmt_hsnow=0.
j = J_jmt
do i=1,imtm2
lcheck = (icy(i,j).or.icy(i+1,j).or.icy(i+1,j-1)
& .or.icy(i+2,j))
area=radius**2*dlambda(i+1)*(1.-sin(phi(j-1)))
facey = radius*cos(phi(j-1))*dlambda(i+1)
vm = vice(i,j-1)
if (vm.ge.0.) qiny = vm*hice_init(i+1,j-1)*hmask(i+1,j-1)
if (vm.lt.0.) qiny = vm*hice_init(i+1,j)*hmask(i+1,j)
if (vm.ge.0.) qay = vm*aice_init(i+1,j-1)*hmask(i+1,j-1)
if (vm.lt.0.) qay = vm*aice_init(i+1,j)*hmask(i+1,j)
if (vm.ge.0.) qsy = qay*hsno_init(i+1,j-1)
if (vm.lt.0.) qsy = qay*hsno_init(i+1,j)
total_jmt_hice=total_jmt_hice+hice_init(i+1,j)
& +qiny*facey*dtts/area
total_jmt_aice=total_jmt_aice+aice_init(i+1,j)
& +qay*facey*dtts/area
if (aice(i+1,j).gt.0.0) then
total_jmt_hsnow=total_jmt_hsnow+hsno_init(i+1,j)
& +(qsy*facey*dtts/area)/aice(i+1,j)
endif
enddo
write(6,*)'area:',area
imtm2r=1./imtm2
total_jmt_hice=total_jmt_hice*imtm2r
total_jmt_aice=total_jmt_aice*imtm2r
total_jmt_hsnow=total_jmt_hsnow*imtm2r
do i=1,imt
hice(i,j)=total_jmt_hice
aice(i,j)=total_jmt_aice
hsnow(i,j)=total_jmt_hsnow
enddo
endif
*DECLARE ICEFLOE
*I ORH3F402.178
if (L_ICYNPOL) then
c SET HEAT FLUXES TO ZERO ON POLAR ISLAND TILL MEANINGFUL
c ONES AER AVAILABLE
ATMSFLUX(i,j_jmt)=0.
SNOWRATE(i,j_jmt)=0.
SUBLIM(i,j_jmt)=0.
CARYHEAT(i,j_jmt)=0.
TOPMELT(i,j_jmt)=0.
BOTMELT(i,j_jmt)=0.
OCEANFLX(i,j_jmt)=0.
HEATFLUX(i,j_jmt)=0.
TOPMELTZ(i,j_jmt)=0.
BOTMELTZ(i,j_jmt)=0.
SUBLIMZ(i,j_jmt)=0.
ATMSFLUX(i,jmt)=0.
SNOWRATE(i,jmt)=0.
SUBLIM(i,jmt)=0.
CARYHEAT(i,jmt)=0.
TOPMELT(i,jmt)=0.
BOTMELT(i,jmt)=0.
OCEANFLX(i,jmt)=0.
HEATFLUX(i,jmt)=0.
TOPMELTZ(i,jmt)=0.
BOTMELTZ(i,jmt)=0.
SUBLIMZ(i,jmt)=0.
endif
*I OCS0F400.13
*IF DEF,MPP
IF (JFIN.GE.JMT_GLOBAL.AND.JST.LE.JMT_GLOBAL) THEN
*ENDIF
if (L_ICYNPOL.AND.(JFIN.EQ.JMT_GLOBAL)) then
do i=1,imt
delh(i,j_jmt)=0.
caryheat(i,j_jmt)=0.
enddo
endif
*IF DEF,MPP
endif
*ENDIF
*DECLARE ICEDRIFT
*I ORH6F402.116
& ,L_ICYNPOL
*I OJG1F405.56
logical icy_old(imt,jmt) ! initial icy
*I ICEDRIFT.349
icy_old(i,j) = (aice_old(i,j).gt.zero)
*D OJG1F405.61
IF (ICY(I,J).or.ICY_OLD(I,J)) THEN
*I OJG1F405.62
ELSE
delhi_dyn(I,J) = 0.
*D OJG1F405.70
IF (ICY(I,J).or.ICY_OLD(I,J)) THEN
*I OJG1F405.71
ELSE
ddt_aice_dyn(i,j)=0.
*D OJG1F405.79
IF (ICY(I,J).or.ICY_OLD(I,J)) THEN
*I OJG1F405.80
ELSE
ddt_hice_dyn(i,j)=0.
*D OJG1F405.88
IF (ICY(I,J).or.ICY_OLD(I,J)) THEN
*I OJG1F405.90
ELSE
delhs_dyn(i,j)=0.
*D OJG1F405.98
IF (ICY(I,J).or.ICY_OLD(I,J)) THEN
*I OJG1F405.100
ELSE
ddt_snow_dyn(i,j)=0.
*/END OF ADVECTION OVER THE POLAR ISLAND SECTION
*DECLARE ICE_CTL
*I JT161193.153
C
real uice_fd(IMT,JMTM1)
real vice_fd(IMT,JMTM1)
real uice_evp_fd(IMT,JMTM1) ! EVP-FD, U
real vice_evp_fd(IMT,JMTM1) ! EVP-FD, V
C
*I OJC3F400.4
&,sig11ne
&,sig11se
&,sig11sw
&,sig11nw
&,sig12ne
&,sig12se
&,sig12sw
&,sig12nw
&,sig22ne
&,sig22se
&,sig22sw
&,sig22nw
*I JT161193.121
&,sstiltx,sstilty
*D ODC1F405.254,ODC1F405.255
&,WSX_ICE(imt_idr,jmtm1_idr) ! IN ZONAL COMPONENT OF WIND STRESS.
&,WSY_ICE(imt_idr,jmtm1_idr) ! IN MERIDIONAL COMPONENT OF " "
*I ODC1F405.257
REAL
& sig11ne(IMT_EVP,JMT_EVP) !Ice internal stresses
&,sig11se(IMT_EVP,JMT_EVP) !for EVP ice dynamics
&,sig11sw(IMT_EVP,JMT_EVP)
&,sig11nw(IMT_EVP,JMT_EVP)
&,sig12ne(IMT_EVP,JMT_EVP)
&,sig12se(IMT_EVP,JMT_EVP)
&,sig12sw(IMT_EVP,JMT_EVP)
&,sig12nw(IMT_EVP,JMT_EVP)
&,sig22ne(IMT_EVP,JMT_EVP)
&,sig22se(IMT_EVP,JMT_EVP)
&,sig22sw(IMT_EVP,JMT_EVP)
&,sig22nw(IMT_EVP,JMT_EVP)
*I ODC1F405.267
&,sstiltx(IMT_EVP,JMTM1_EVP),sstilty(IMT_EVP,JMTM1_EVP)
*D ODC1F405.287,ORH1F305.5717
*D ODC1F405.289,ODC1F405.303
if (L_ICEEVP) THEN
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
&,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_ice,wsy_ice,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
*D JT161193.208,ODC1F405.324
IF (L_OTIMER) CALL TIMER('ICE_EVP',4)
C
ENDIF !L_ICEEVP = true
*/
*/ modify icedrift call for ICESIMPLE
*/
*I ODC1F405.331
C
CALL ICEDRIFT(
*CALL ARGOINDX
C & AICE,HICE,HSNOW,uice_local,vice_loc
C WMC make sure we're using the right ice velocity...
& AICE,HICE,HSNOW,uice ,vice ,ICY,NEWICE
&,OCEAN,IMT,IMTM1,IMTM2,JMT,JMTM1,KM,FKMP
&,PHI,PHIT,AICEMIN,AMX,HICEMIN,ah_ice,DPHI,DLAMBDA
+,DTTS
+,CS,CSTR,DYUR,DYTR,DXU2R,DXT4R
&,RADIUS_SI
+,STASHWORK(SI(202,32,im_index)),STASHWORK(SI(203,32,im_index))
+,STASHWORK(SI(209,32,im_index)),STASHWORK(SI(210,32,im_index))
+,STASHWORK(SI(201,32,im_index)),STASHWORK(SI(204,32,im_index))
+,SF(202,32),SF(203,32),SF(209,32),SF(210,32)
+,SF(201,32),SF(204,32)
&,stashwork(si(223,32,im_index)),stashwork(si(224,32,im_index))
&,stashwork(si(225,32,im_index)),stashwork(si(226,32,im_index))
&,sf(223,32),sf(224,32),sf(225,32),sf(226,32)
+ )
C
*/
*/ modify icedrift call for EVP
*/
*D ODC1F405.333
IF (L_ICEEVP) THEN
*D ODC1F405.337
& AICE,HICE,HSNOW,UICE,VICE,ICY,NEWICE
*/
*/*I ICE_CTL.287
*DECK ICEEVP
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_ice,wsy_ice,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
& )
implicit none
*CALL CSUBMODL
*CALL PARPARM
*CALL TYPSIZE
*CALL TYPSTS
*CALL TYPIEVP
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)
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_ice,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
c CALL swap_1d(wsy_ice,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
& )
do k=1,ndte ! subcycling
call stress( ! stress tensor equation (sigma)
*CALL ARGIEVP
&,coslamda,sinlamda
&,JFIN,JMT_GLOBAL
& )
call stepu( ! total surface stress, momentum eqn
*CALL ARGIEVP
&,uice_init,vice_init
&,dsig11dx
&,dsig12dy
&,dsig21dx
&,dsig22dy
&,JFIN,JMT_GLOBAL
&,coslamda,sinlamda
& )
enddo
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)
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
& xymin_buffer !store for calculation of xymin
&,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
cc calculate dte,dter,ecc2m,ecc2p,
cc ocean_mask,t_dx,t_dy,t_dxr,t_dyr,HT_N,HT_E,xymin
cc Subcycle timestep variables
c dte=dtts/ndte
dter=1.0/dte
ndte=dtts*dter
tdamp2=2.*eyc*dtts
C
IF (first_time_called) THEN
write(6,*)'Namelist parameters from ICEEVPNL'
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,*)'dte=',dte
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
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
do j=1,jmt
do i=1,imt
tmass(i,j) = rhoice*hice(i,j) ! k/m^2
& +rhosnow*aice(i,j)*hsnow(i,j)
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
if (L_ICESSTILT) then
do j=J_1,J_JMTM1
do i=1,imt
wsx_ice(i,j) = ocean_uv_mask(i,j)*
& (wsx_ice(i,j) - g*umass(i,j)*sstiltx(i,j))
wsy_ice(i,j) = ocean_uv_mask(i,j)*
& (wsy_ice(i,j) - g*umass(i,j)*sstilty(i,j))
enddo
enddo
else
do j=J_1,J_JMTM1
do i=1,imt
wsx_ice(i,j) = ocean_uv_mask(i,j)*
& wsx_ice(i,j)
wsy_ice(i,j) = ocean_uv_mask(i,j)*
& wsy_ice(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
& )
c.. Calculates the internal stress components, sigma_ij, in the four
c.. triangles of each cell.
implicit none
*CALL TYPIEVP
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.. 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
c if((JFIN .ge. JMT_GLOBAL).and.L_ICYNPOL) then
cc FIND A SOLUTION ON THE POLAR ISLAND
cc CALCULATE STRAIN RATE ON LOCALLY CARTESIAN COORDS
cc (1: perpendicular to 2, and horizontal,2:The Greenwich meridian)
c xi11pol=0.
c xi12pol=0.
c xi22pol=0.
c j=j_jmt
c do i=2,imtm1
c xi11pol=xi11pol+sinlamda(i)*(uice(i,j-1)*coslamda(i)-
c & vice(i,j-1)*sinlamda(i))
c & *0.5*t_dyr(i,j)
c xi22pol=xi22pol-coslamda(i)*(uice(i,j-1)*sinlamda(i)+
c & vice(i,j-1)*coslamda(i))
c & *0.5*t_dyr(i,j)
c xi12pol=xi12pol+uice(i,j-1)*0.5*t_dyr(i,j)
c enddo
c deltapol = sqrt( (xi11pol**2+xi22pol**2)*eccp
c & + xi12pol**2*ecci + xi11pol*xi22pol*eccm)
c prsspol = pstar*hice(1,j)*exp(-cstar*(1.-aice(1,j)))
c prsdpol=prsspol/max(rcon,deltapol)
c prs2pol=0.5*prsdpol*deltapol
cc For polar point solution all ne/nw/se/sw triangles are equated
cc and each point i=1 to imt is equal but rotated to be locally
cc Cartesian with i=1 as a basis. So will use sig**sw(i=1,j_jmt)
cc as history in time derivatives
cC CHECK THIS IS CONSISTENT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
c c4pol=Tdter*sig11sw(1,j_jmt)+prsdpol*xi11pol-prs2pol
c c5pol=Tdter*sig22sw(1,j_jmt)+prsdpol*xi22pol-prs2pol
c sig11pol=acoef*c5pol+bcoef*c4pol
c sig22pol=acoef*c4pol+bcoef*c5pol
c sig12pol=(Tdter*sig12sw(1,j_jmt)+prsdpol*xi12pol)*ccoef
c do i=1,imt
c sig11ne(i,j_jmt)=coslamda(i)*sig11pol+sinlamda(i)*sig22pol
c sig22ne(i,j_jmt)=-sinlamda(i)*sig11pol+coslamda(i)*sig22pol
c sig12ne(i,j_jmt)=sig12pol
c sig11nw(i,j_jmt)=sig11ne(i,j_jmt)
c sig11se(i,j_jmt)=sig11ne(i,j_jmt)
c sig11sw(i,j_jmt)=sig11ne(i,j_jmt)
c sig12nw(i,j_jmt)=sig12ne(i,j_jmt)
c sig12se(i,j_jmt)=sig12ne(i,j_jmt)
c sig12sw(i,j_jmt)=sig12ne(i,j_jmt)
c sig22nw(i,j_jmt)=sig22ne(i,j_jmt)
c sig22se(i,j_jmt)=sig22ne(i,j_jmt)
c sig22sw(i,j_jmt)=sig22ne(i,j_jmt)
c enddo
c write(6,*)'prsspol: ',prsspol
c write(6,*)'xi11pol: ',xi11pol
c write(6,*)'xi12pol: ',xi12pol
c write(6,*)'xi22pol: ',xi22pol
c write(6,*)'sig11pol: ',sig11pol
c write(6,*)'sig12pol: ',sig12pol
c write(6,*)'sig22pol: ',sig22pol
c endif
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_ice(i,j_jmt)=wsx_ice(i,j_jmt-1)
wsy_ice(i,j_jmt)=wsy_ice(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
taux = wsx_ice(i,j) + isx(i,j)
c & +tiltx(i,j)
tauy = wsy_ice(i,j) + isy(i,j)
c & +tilty(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
*/---------------------------------------------------------------------
@
1.1
log
@Initial revision
@
text
@d12 1
d706 3
a708 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
d712 3
a714 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
d718 3
a720 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
d724 3
a726 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
d731 3
a733 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
d737 3
a739 1
CALL COPYODIAGL(IMT,JMT,nlev,.TRUE.,0.
@