From 02fd787f7ab691e753561326719b7d886f5ff1b2 Mon Sep 17 00:00:00 2001 From: wangqian2284 <369402617@qq.com> Date: Tue, 26 Mar 2024 02:17:20 +0800 Subject: [PATCH] Update mice_mevp.F90 add ics==1 in mice_mevp for xy coordinate --- src/Multi_ice/mice_mevp.F90 | 92 +++++++++++++++++++++++-------------- 1 file changed, 58 insertions(+), 34 deletions(-) diff --git a/src/Multi_ice/mice_mevp.F90 b/src/Multi_ice/mice_mevp.F90 index 7fb07a332..4c59f8b7d 100644 --- a/src/Multi_ice/mice_mevp.F90 +++ b/src/Multi_ice/mice_mevp.F90 @@ -1,7 +1,7 @@ ! Modified EVP dynamics of ice module (Bouillon 2013) subroutine ice_mevp use schism_glbl,only: rkind,time_stamp,eta2,np,npa,ne,nea,dldxy,elnode,i34,cori, & - &grav,isbnd,nne,indel,area,iself,time_stamp,rnday,fdb,lfdb,xnd,ynd,iplg,ielg, & + &grav,isbnd,nne,indel,area,iself,time_stamp,rnday,fdb,lfdb,xnd,ynd,iplg,ielg,ics, & &elside,mnei,rho0,idry,errmsg,pframe,eframe,indnd,nnp,omega_e,xlon,ylat,idry_e,pr use schism_msgp, only: myrank,nproc,parallel_abort,exchange_p2d use mice_module @@ -16,13 +16,13 @@ subroutine ice_mevp real(rkind) :: tmp1,tmp2,pp0,delta,rr1,rr2,rr3,sig1,sig2,x10,x20,y10,y20,rl10, & &rl20,sintheta,bb1,bb2,h_ice_el,a_ice_el,h_snow_el,dsig_1,dsig_2,mass, & &cori_nd,umod,gam1,rx,ry,rxp,ryp,eps11,eps12,eps22, & - &zeta,delta_nd,sum1,sum2,dt_by_mass,tmp3,tmp4 + &zeta,delta_nd,sum1,sum2,dt_by_mass,tmp3,tmp4,beta0,mevp_alpha3,mevp_alpha4 integer :: iball(mnei) real(rkind) :: swild(2,3),deta(2,nea), & &swild2(nea),alow(4),bdia(4),rrhs(3,4),U_ice_0(npa),V_ice_0(npa),utmp(3),vtmp(3), & &a_ice0_0(npa),a_icen0(npa,ncat),v_icen0(npa,ncat),a_icen_elem(ncat),v_icen_elem(ncat),swild1(3),& - &deta_pice(2,nea),p_ice(3) + &deta_pice(2,nea),p_ice(3),al_beta(nea) logical :: & calc_strair @@ -34,6 +34,8 @@ subroutine ice_mevp a_ice0_0(:) = 0 a_icen0(:,:) = 0 v_icen0(:,:) = 0 + mevp_alpha3 = 200 + mevp_alpha4 = 0.02 call icepack_to_schism (nx_in=npa, & aice_out=a_ice0, & @@ -42,6 +44,11 @@ subroutine ice_mevp aice0_out=a_ice0_0, & aicen_out=a_icen0, & vicen_out=v_icen0 ) + + + do i=1,nea + al_beta(i)=mevp_alpha3/tanh(mevp_alpha4*area(i)/dt_dyn) + enddo !i=1,nea !do i = 1,nea ! swild1(1) = sum(a_ice0(elnode(1:i34(i),i)))/i34(i) ! swild1(2) = sum(m_ice0(elnode(1:i34(i),i)))/i34(i) @@ -93,12 +100,13 @@ subroutine ice_mevp tmp3 = U_ice(elnode(j,i)) tmp4 = V_ice(elnode(j,i)) - !utmp(j) = tmp3*tmp1+tmp2*tmp4 - !vtmp(j) = tmp1*tmp4+tmp3*tmp2 - !utmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,1,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,1,i)) - !vtmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,2,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,2,i)) - utmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,1,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,1,i)) - vtmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,2,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,2,i)) + if(ics==1) then + utmp(j) = tmp3 + vtmp(j) = tmp4 + else + utmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,1,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,1,i)) + vtmp(j) = tmp3*dot_product(pframe(1:3,1,elnode(j,i)),eframe(1:3,2,i))+tmp4*dot_product(pframe(1:3,2,elnode(j,i)),eframe(1:3,2,i)) + endif enddo eps11=dot_product(utmp,dldxy(1:i34(i),1,i)) !epsilon11=du_dx eps22=dot_product(vtmp,dldxy(1:i34(i),2,i)) !epsilon22=dv_dy @@ -128,8 +136,11 @@ subroutine ice_mevp sig1=sigma11(i)+sigma22(i) !from previous iteration sig2=sigma11(i)-sigma22(i) - sig1=sig1+(rr1-sig1)/mevp_alpha1 - sig2=sig2+(rr2-sig2)/mevp_alpha1 + !sig1=sig1+(rr1-sig1)/mevp_alpha1 + !sig2=sig2+(rr2-sig2)/mevp_alpha1 + + sig1=sig1+(rr1-sig1)/al_beta(i) + sig2=sig2+(rr2-sig2)/al_beta(i) sigma12(i)=sigma12(i)+(rr3-sigma12(i))/mevp_alpha1 sigma11(i)=0.5*(sig1+sig2) @@ -170,7 +181,11 @@ subroutine ice_mevp else if(maxval(idry(elnode(1:i34(i),i)))/=0) then !dry deta(1:2,i)=0 - deta_pice(1:2,i) = 0 + p_ice=(rhoice*m_ice0(elnode(1:i34(i),i))+rhosno*m_snow0(elnode(1:i34(i),i))+pr(elnode(1:i34(i),i))/grav)/rho0 + deta_pice(1,i)=dot_product(p_ice(1:3),dldxy(1:i34(i),1,i)) + deta_pice(2,i)=dot_product(p_ice(1:3),dldxy(1:i34(i),2,i)) + !deta_pice(1,i) = 0 !dot_product((p_ice(1:3)),dldxy(1:i34(i),1,i)) + !deta_pice(2,i) = 0 !dot_product((p_ice(1:3)),dldxy(1:i34(i),2,i)) else !wet p_ice=(rhoice*m_ice0(elnode(1:i34(i),i))+rhosno*m_snow0(elnode(1:i34(i),i))+pr(elnode(1:i34(i),i))/grav)/rho0 deta_pice(1,i)=dot_product((eta2(elnode(1:i34(i),i))+p_ice(1:3)),dldxy(1:i34(i),1,i)) @@ -195,10 +210,14 @@ subroutine ice_mevp cycle endif - if(idry(i)/=0) then !dry - U_ice(i)=0; V_ice(i)=0 - cycle - endif + !if(idry(i)/=0) then !dry + !if(maxval(idry(indnd(1:nnp(i),i)))>0) then + ! U_ice(i)=0; V_ice(i)=0 + ! cycle + !endif + + beta0=dot_product(weit_elem2node(1:nne(i),i),al_beta(iball(1:nne(i)))) + !Not bnd node; has ice mass=(rhoice*ice_tr0(1,i)+rhosno*ice_tr0(3,i)) !>0 mass=max(mass,9.d0*ice_tr0(2,i)) !limit m/a>=9 @@ -218,16 +237,16 @@ subroutine ice_mevp !gam1=ice_tr0(2,i)*dt_by_mass*cdwat*rhowat*umod if (calc_strair) then - rx=mevp_alpha2*U_ice(i)+U_ice_0(i)+gam1*(u_ocean(i)*cos_io-v_ocean(i)*sin_io)+ & + rx=beta0*U_ice(i)+U_ice_0(i)+gam1*(u_ocean(i)*cos_io-v_ocean(i)*sin_io)+ & &dt_by_mass*stress_atmice_x(i) - ry=mevp_alpha2*V_ice(i)+V_ice_0(i)+gam1*(u_ocean(i)*sin_io+v_ocean(i)*cos_io)+ & + ry=beta0*V_ice(i)+V_ice_0(i)+gam1*(u_ocean(i)*sin_io+v_ocean(i)*cos_io)+ & &dt_by_mass*stress_atmice_y(i) else - rx=mevp_alpha2*U_ice(i)+U_ice_0(i)+gam1*(u_ocean(i)*cos_io-v_ocean(i)*sin_io)+ & + rx=beta0*U_ice(i)+U_ice_0(i)+gam1*(u_ocean(i)*cos_io-v_ocean(i)*sin_io)+ & &dt_by_mass*ice_tr0(2,i)*stress_atmice_x(i) - ry=mevp_alpha2*V_ice(i)+V_ice_0(i)+gam1*(u_ocean(i)*sin_io+v_ocean(i)*cos_io)+ & + ry=beta0*V_ice(i)+V_ice_0(i)+gam1*(u_ocean(i)*sin_io+v_ocean(i)*cos_io)+ & &dt_by_mass*ice_tr0(2,i)*stress_atmice_y(i) - endif + endif !Pressure gradient sum1=0; sum2=0 @@ -236,13 +255,15 @@ subroutine ice_mevp h_ice_el=sum(ice_tr0(1,elnode(1:3,ie)))/3.0 h_snow_el=sum(ice_tr0(3,elnode(1:3,ie)))/3.0 tmp2 = rhoice*h_ice_el+rhosno*h_snow_el !mass @elem - tmp2 = mass - !sum1=sum1+tmp2*deta(1,ie)*area(ie)/3 - !sum2=sum2+tmp2*deta(2,ie)*area(ie)/3 - - sum1=sum1+tmp2*deta_pice(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,1,i))+tmp2*deta_pice(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,1,i)) - sum2=sum2+tmp2*deta_pice(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,2,i))+tmp2*deta_pice(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,2,i)) + !tmp2 = mass + if(ics==1) then + sum1=sum1+tmp2*deta(1,ie)*area(ie)/3 + sum2=sum2+tmp2*deta(2,ie)*area(ie)/3 + else + sum1=sum1+tmp2*deta_pice(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,1,i))+tmp2*deta_pice(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,1,i)) + sum2=sum2+tmp2*deta_pice(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,2,i))+tmp2*deta_pice(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,2,i)) + endif !sum1=sum1+tmp2*deta(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,1,i))+tmp2*deta(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,1,i)) !sum2=sum2+tmp2*deta(1,ie)*area(ie)/3*dot_product(eframe(1:3,1,ie),pframe(1:3,2,i))+tmp2*deta(2,ie)*area(ie)/3*dot_product(eframe(1:3,2,ie),pframe(1:3,2,i)) enddo !j @@ -254,17 +275,20 @@ subroutine ice_mevp do j=1,nne(i) ie=indel(j,i) id=iself(j,i) - !sum1=sum1+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie)) - !sum2=sum2+area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie)) - sum1=sum1+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie))*dot_product(eframe(1:3,1,ie),pframe(1:3,1,i))+& - &area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie))*dot_product(eframe(1:3,2,ie),pframe(1:3,1,i)) - sum2=sum2+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie))*dot_product(eframe(1:3,1,ie),pframe(1:3,2,i))+& - &area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie))*dot_product(eframe(1:3,2,ie),pframe(1:3,2,i)) + if(ics==1) then + sum1=sum1+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie)) + sum2=sum2+area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie)) + else + sum1=sum1+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie))*dot_product(eframe(1:3,1,ie),pframe(1:3,1,i))+& + &area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie))*dot_product(eframe(1:3,2,ie),pframe(1:3,1,i)) + sum2=sum2+area(ie)*(dldxy(id,1,ie)*sigma11(ie)+dldxy(id,2,ie)*sigma12(ie))*dot_product(eframe(1:3,1,ie),pframe(1:3,2,i))+& + &area(ie)*(dldxy(id,1,ie)*sigma12(ie)+dldxy(id,2,ie)*sigma22(ie))*dot_product(eframe(1:3,2,ie),pframe(1:3,2,i)) + endif enddo !j rx=rx-dt_by_mass/area_median(i)*sum1 ry=ry-dt_by_mass/area_median(i)*sum2 - tmp1=1+mevp_alpha2+gam1*cos_io+Tbu(i)/(sqrt(U_ice_0(i)**2+V_ice_0(i)**2)+5e-5) + tmp1=1+beta0+gam1*cos_io+Tbu(i)/(sqrt(U_ice_0(i)**2+V_ice_0(i)**2)+5e-5) tmp2=dt_dyn*cori_nd+gam1*sin_io delta=tmp1*tmp1+tmp2*tmp2 if(delta<=0) call parallel_abort('ice_mevp: delta<=0')