module driver use const use comm use var use io use util use solver use physics use filter use shr_flux use hmix_GM_submeso_share use hmix_GM implicit real*8 (a-h,o-z) contains subroutine time_management_timcom(itf) use time_management, only: iyear,imonth !Derive several time indicies for loop control, historical avg. output !and seasonal data interpolation. Here we use synergizing time:syng_day to simplify the calulation and !find the current month. !=================================================================================== !derived variables listed below, !days: modeling days (#iteration/ iteration per day) !syng_days: synergy days !nyr: modeling years !nmon: current month indice !nopt: current monthly days number !wopt: weighting for historical monthly avg. output !nld,new: month indicies for seasonal data interpolation !fold,fnew: interpolation weighting for seasonal data integer, intent(in) :: itf !interation loop !set modeling days and years days=itf/daodt n=int(days)/365 nyr=n+1 !find simplified days for later loop control. syng_days=mod(days,365.d0) !find current month dtmp=syng_days nld=0 if(dtmp.eq.0) nld=12 do while(dtmp.gt.0) nld=nld+1 dtmp=dtmp-dmonth(nld) enddo if(nld.eq.0) then dtmp=-1.d0 nld=1 new=2 elseif (nld.eq.12) then new=1 else new=nld+1 endif fold=-1.d0*dtmp/dmonth(nld) fnew=1.d0-fold nmon=nld !hist time avg. variables if (lrstrt==0.and.iyear==1.and.imonth==1) then nopt=dmonth(nmon)-1 else nopt=dmonth(nmon) endif wopt=1./dble(nopt*itfday) end subroutine time_management_timcom subroutine fsglo !use time_management, only:check_time_flag,nsteps_run !use ice, only:ice_cpl_flag character fn1*7 integer tt0,tt1,tt2 dimension rhomin(k1) dimension scr1(i0,j0,k1) real*8, dimension(i0,j0,0:k0,2) :: vdc real*8, dimension(i0,j0) :: works, workt real*8, dimension(k1) :: tzt, tzb,lttt,lbbb,txxx,www,ttt,GGG,GGGG real*8, dimension(i2,j2) :: tz1,t2sst !----check salinity !if (.false.) then totals=sum(s2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) dtmp=totals call mpi_allreduce(dtmp,totals,1,mpi_real8,mpi_sum,m_cart,ierr) totals_sfc=sum(s2(2:i1,2:j1,1:5)*volume(2:i1,2:j1,1:5)) dtmp=totals_sfc call mpi_allreduce(dtmp,totals_sfc,1,mpi_real8,mpi_sum,m_cart,ierr) totalt=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) dtmp=totalt call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) totalt_sfc=sum(t2(2:i1,2:j1,1:5)*volume(2:i1,2:j1,1:5)) dtmp=totalt_sfc call mpi_allreduce(dtmp,totalt_sfc,1,mpi_real8,mpi_sum,m_cart,ierr) ! if (myid==0 .and. mod(ITF,ITFDAY)==0) write(20,*)'meanS',TOTALS/TOTALVOL ! if (myid==0 .and. mod(ITF,ITFDAY)==0) write(21,*)'meanS',TOTALS_sfc/TOTALVOL_sfc ! if (myid==0 .and. mod(ITF,ITFDAY)==0) write(23,*)'meanT',TOTALT_sfc/TOTALVOL_sfc if (mod(ITF,ITFDAY)==1) then !WRITE(*,*)'CHECK--TIMCOM--FORCING',ITF, ITFDAY if (myid==0) write(*,*) 'calculating turbulent fluxes and radiative flux' timer(14)=mpi_wtime() call flux_atmocn timer(14)=mpi_wtime()-timer(14) end if !if (myid .eq. 0) write(*,*) 'myid0', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 1) write(*,*) 'myid1', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 2) write(*,*) 'myid2', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 32) write(*,*) 'myid1', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 33) write(*,*) 'myid1', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 34) write(*,*) 'myid1', xdeg(1), xdeg(10), ydeg(1), ydeg(10) !if (myid .eq. 64) write(*,*) 'myid1', xdeg(1), xdeg(10), ydeg(1), ydeg(10) ! totalt=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt ! call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT01',TOTALT/TOTALVOL if (lkppmix) then ! if (mod(itf,5).eq.1) then ! timer(11)=mpi_wtime() call set_STF call kppglo(vdc) ! timer(11)=mpi_wtime()-timer(11) ! endif else if (mod(itf,5).eq.1) then timer(11)=mpi_wtime() call pp82glonew timer(11)=mpi_wtime()-timer(11) endif endif ! totalt=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt ! call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT02',TOTALT/TOTALVOL call solve_hydrostatic_equation ! modified pacanowski and philander vertical mixing coefficients ! subracting layer min rho to reduce roundoff error cannot be ! done before pp82glo or loop ! if (mod(itf,4).eq.1) call pp82glo ! if (mod(itf,4).eq.1) call pp82glonew if (smagmix) then ! smagorinsky mixing, not useful or needed due to model robustness if (mod(itf,5).eq.1) call smagglo endif ! totalt1=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt1 ! call mpi_allreduce(dtmp,totalt1,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT1',TOTALT1/TOTALVOL call mpi_barrier(m_cart,ierr) timer(6)=mpi_wtime() call mpi_exch_3d_r8(u2) call mpi_exch_3d_r8(t2) call mpi_exch_3d_r8(s2) call mpi_exch_3d_r8(u1) call mpi_exch_3d_r8(t1) call mpi_exch_3d_r8(s1) call mpi_exch_3d_r8(ulf) call mpi_exch_3d_r8(tlf) call mpi_exch_3d_r8(slf) call mpi_exch_3d_r8(p) call mpi_exch_3d_r8(c2) call mpi_exch_3d_r8(c1) call mpi_exch_3d_r8(clf) if (lopen==2) then call mpi_exch_3d_r8_vsymtry(v2) call mpi_exch_3d_r8_vsymtry(v1) call mpi_exch_3d_r8_vsymtry(vlf) else call mpi_exch_3d_r8(v2) call mpi_exch_3d_r8(v1) call mpi_exch_3d_r8(vlf) endif timer(6)=mpi_wtime()-timer(6) lb=1 lt=2 !WRITE(*,*)'lfsrf=',lfsrf if (lfsrf .eq. 0) then ! write(42,*) 'flag1' !do j=1,j2 ! do i=1,i2 ! uz(i,j,1)=0. ! vz(i,j,1)=0. ! sz(i,j,1)=0. ! tz(i,j,1)=0. ! cfcz(i,j,1)=0. ! enddo !enddo uz(1:i2, 1:j2, 1) = 0. vz(1:i2, 1:j2, 1) = 0. sz(1:i2, 1:j2, 1) = 0. tz(1:i2, 1:j2, 1) = 0. cfcz(1:i2, 1:j2, 1) = 0. else !write(42,*) 'flag2' !do j=1,j2 ! do i=1,i2 ! uz(i,j,1)=w(i+1,j+1,1)*u2(i+1,j+1,1) ! vz(i,j,1)=w(i+1,j+1,1)*v2(i+1,j+1,1) ! sz(i,j,1)=w(i+1,j+1,1)*s2(i+1,j+1,1) ! tz(i,j,1)=w(i+1,j+1,1)*t2(i+1,j+1,1) ! cfcz(i,j,1)=w(i+1,j+1,1)*c2(i+1,j+1,1) ! enddo !enddo uz(1:i2, 1:j2, 1) = w(2:i2+1, 2:j2+1, 1) * u2(2:i2+1, 2:j2+1,1) vz(1:i2, 1:j2, 1) = w(2:i2+1, 2:j2+1, 1) * v2(2:i2+1, 2:j2+1,1) sz(1:i2, 1:j2, 1) = w(2:i2+1, 2:j2+1, 1) * s2(2:i2+1, 2:j2+1,1) tz(1:i2, 1:j2, 1) = w(2:i2+1, 2:j2+1, 1) * t2(2:i2+1, 2:j2+1,1) !2021/9/13/ 2pm test ! t2sst(1:i2,1:j2) =t2(2:i2+1, 2:j2+1, 1) ! tz1(1:i2,1:j2) = tz(1:i2,1:j2,1) !20210915 cfcz(1:i2, 1:j2, 1) = w(2:i2+1, 2:j2+1, 1) * c2(2:i2+1, 2:j2+1, 1) endif ! totalt111=sum(t1(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt111 ! call mpi_allreduce(dtmp,totalt111,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT11*',TOTALT111/TOTALVOL ! ================================================================= !calculate all internal fluxes and substitute into conservation laws ! ================================================================= ! level-by-level moving window reduces storage. ! note: all boundary fluxes are masked out in this loop, ! but are included later. this is the main computation loop. do k=1,k1 l=k+1 ! variable diffusivity arrays are used for near-shore laplace filter !do j=1,j0 ! do i=1,i1 ! dhx(i,j)=dmx(i,j,k) ! enddo !enddo dhx(1:i1, 1:j0) = dmx(1:i1, 1:j0, k) !dhx(1:i1, 1:j0) = .1d0*dmx(1:i1, 1:j0, k) !do j=1,j1 ! do i=1,i0 ! dhy(i,j)=dmy(i,j,k) ! enddo !enddo dhy(1:i0, 1:j1) = dmy(1:i0, 1:j1, k) !dhy(1:i0, 1:j1) = .1d0*dmy(1:i0, 1:j1, k) !================================================================ ! (Sep 2018) ! Gent McWilliams (GM) parameterization ! for eddy-induced tracer transport and diffusion in ocean !================================================================ if ( hmix_tracer_type_gm > 0.5) then if (k == 1) then call init_meso_mixing call init_gm call tracer_diffs_and_isopyc_slopes(trcr) endif call hdifft_gm (k, trcr, vdc, HOR_DIFF_X, HOR_DIFF_Y, KH_GM_X, KH_GM_Y, GTK0) GGG(k)=sum(GTK0(2:i1,2:j1,1)) ! GGGG(k)=GTK0(2,2,1)! if (k0 --> dval>0 tnud_tmp(i-1,j-1)=in(i,j,1)*t_nudstr*(tkli-t2(i,j,1)) snud_tmp(i-1,j-1)=in(i,j,1)*s_nudstr*(skli-s2(i,j,1)) ! accumulate nudgings from present month for later qavg calc. tnudge(i-1,j-1)=tnudge(i-1,j-1)+tnud_tmp(i-1,j-1) snudge(i-1,j-1)=snudge(i-1,j-1)+snud_tmp(i-1,j-1) !first year use initialized qavg value. if (days <=365.) then tnudge=0.d0 snudge=0.d0 endif ! qavg(watts/cm-cm)= (watts/deg/cc)*heat*dx(j)*dy*dz/(dx(j)*dy)/time ! where qavg=heat/time in deg ! per time step !====================================================================== ! surface restoring: nudge toward surface t climatology ! to emulate eddy damping by atmospheric heat flux. !tnud_tmp(i-1,j-1)=( in(i,j,1)*(tmp2*t2(i,j,1)+tmp1*tkli)+ & ! (1-in(i,j,1))*tkli )-t2(i,j,1) enddo enddo !calculating model determined monthly surface fluxes term at the end of the month, !first year use initialized value if(abs(syng_days-epday(nmon)).gt.1.e-6 .or. days .le. 365. ) go to 575 if(myid==0) write(*,*) 'deriving monthly determined ensemble fluxes' temp=1./(dmonth(nmon)*itfday) nsombo(nmon)= nsombo(nmon)+1 !taus=(30.d0*24.d0*3600.d0) ! 30days !tmp2=z(3)*ens1/taus do j=1,j2 do i=1,i2 ! EVAP IS THE NET SURFACE OUTFLOW (POSITIVE FOR NET POSITIVE E-P) ! WAVG=-(DZ*DS/S)/DT, SELINC: CHANGES OF SALINITY PER TIME STEP !selinc=-temp*(s2(i+1,j+1,1)-sclis(i+1,j+1,1,new)) !wavg(i,j,nld)=(wavg(i,j,nld)-tmp2*selinc/s2(i+1,j+1,1))*IN(i+1,j+1,1) !wtmp(i,j,nld)=wtmp(i,j,nld)+wavg(i,j,nld)+ & ! tmp2*(s2(i+1,j+1,1)-sclis(i+1,j+1,1,new))/s2(i+1,j+1,1) wtmp(i,j,nld)=wtmp(i,j,nld)+wavg(i,j,nld)+ & temp*(snudge(i,j)+sclis(i+1,j+1,1,new)-s2(i+1,j+1,1)) qtmp(i,j,nld)=qtmp(i,j,nld)+qavg(i,j,nld)+ & temp*(tnudge(i,j)+tclis(i+1,j+1,1,new)-t2(i+1,j+1,1)) qavg(i,j,nld)=in(i+1,j+1,1)*qtmp(i,j,nld)/nsombo(nmon) wavg(i,j,nld)=in(i+1,j+1,1)*wtmp(i,j,nld)/nsombo(nmon) ! Initialize accumulated nudgings for the next month tnudge(i,j)=0.d0 snudge(i,j)=0.d0 enddo enddo if (days<=365.)then ! ====================================================================== ! first year only (gives reasonable advection estimate for next month; ! also, avoids accumulation of advection-induced drift do j=1,j2 do i=1,i2 qavg(i,j,new)=qavg(i,j,nld) wavg(i,j,new)=wavg(i,j,nld) enddo enddo endif !Finally, update the restoring and nudging fluxes. !heat flux. 575 if(myid==0) write(*,*) 'update restoring fluxes' !compute epcomp for global fresh water conservation !dSall=runoff,evap,prec,salt,wavg,snudging !dA(1)=total ocean surface area !epcomp=sum(dSall*dA(1))/sum(dA(1)) call fw_conserv if(.false.)then do j=2,j1 ccd=dx(j)*dy(j) do i=2,i1 !temperature temp=in(i,j,1)*qavg(i-1,j-1,nld) t2(i,j,1)=t2(i,j,1)+temp+tnud_tmp(i-1,j-1) t1(i,j,1)=t1(i,j,1)+temp+tnud_tmp(i-1,j-1) tlf(i,j,1)=tlf(i,j,1)+temp+tnud_tmp(i-1,j-1) !Salinity temp=in(i,j,1)*(wavg(i-1,j-1,nld)-epcomp) s2(i,j,1)=s2(i,j,1)+temp+snud_tmp(i-1,j-1) s1(i,j,1)=s1(i,j,1)+temp+snud_tmp(i-1,j-1) slf(i,j,1)=slf(i,j,1)+temp+snud_tmp(i-1,j-1) enddo enddo endif ! ============= ! bottom stress ! ============= tmp=dt*drag do j=2,j1 do i=2,i1 k=kb(i,j) ! don't exceed explicit limit drg=min(.2,in(i,j,k)*tmp*odz(k)*sqrt(u1(i,j,k)**2+v1(i,j,k)**2)) u2(i,j,k)=u2(i,j,k)-drg*u1(i,j,k) v2(i,j,k)=v2(i,j,k)-drg*v1(i,j,k) enddo enddo ! ==================== ! trapezoidal coriolis ! ==================== do k=1,k1 do j=2,j1 do i=2,i1 tmp=.5*(f(j-1)+ulf(i,j,k)*tanphi(j-1))*dt temp=1./(1.+tmp**2) qu=u2(i,j,k)+tmp*v1(i,j,k) qv=v2(i,j,k)-tmp*u1(i,j,k) u2(i,j,k)=temp*(qu+tmp*qv) v2(i,j,k)=temp*(qv-tmp*qu) enddo enddo enddo ! ====================================== ! open boundary conditions at North Pole ! ====================================== timer(7)=mpi_wtime() if (lopen==2) then ! adapt Symmetric boundary conditions at North Pole call mpi_exch_3d_r8(u2) call mpi_exch_3d_r8_vsymtry(v2) else ! nbv update must be after loop 500 (to keep non-divergent bt mode) if (lopen==1 .and. mylat .eq. ngy1) then tmp=.5*dt do k=1,k1 ! no open boundaries in periodic longitudinal direction do i=2,i1 tmpin=in(i,j1,k)*tmp*ody(j1)*csv(j1)*ocs(j1) temp=abs(v(i,j1,k)) temp1=tmpin*(temp+v(i,j1,k)) temp2=tmpin*(temp-v(i,j1,k)) u2(i,j1,k)=u2(i,j1,k)-temp1*ulf(i,j1,k)+temp2*u2(i,j0,k) v2(i,j1,k)=v2(i,j1,k)-temp1*vlf(i,j1,k)+temp2*v2(i,j0,k) s2(i,j1,k)=s2(i,j1,k)-temp1*slf(i,j1,k)+temp2*s2(i,j0,k) t2(i,j1,k)=t2(i,j1,k)-temp1*tlf(i,j1,k)+temp2*t2(i,j0,k) c2(i,j1,k)=c2(i,j1,k)-temp1*clf(i,j1,k)+temp2*c2(i,j0,k) enddo enddo endif call mpi_exch_3d_r8(u2) call mpi_exch_3d_r8(v2) endif timer(7)=mpi_wtime()-timer(7) call mpi_barrier(m_cart,ierr) ! ===================================================== ! interpolate "a" grid quantities to "c" grid locations ! ===================================================== ! 4th-order interpolations ! unnecessary masking has been removed here do k=1,k1 do j=2,j1 !do i=1,i1 ! scr(i,j,k)=6.*(u2(i,j,k)+u2(i+1,j,k)) !enddo scr(1:i1, j, k) = 6. * (u2(1:i1, j, k) + u2(2:i1+1, j, k)) ! skip loop 612 for 2nd order scheme do i=2,i2 scr(i,j,k)=scr(i,j,k)+iu(i,j,k)*(-u2(i-1,j,k)+u2(i,j,k)+u2(i+1,j,k)-u2(i+2,j,k)) enddo !scr(2:i2,j,k) = scr(2:i2,j,k) + iu(2:i2,j,k) * (-u2(1:i2-1,j,k)+u2(2:i2,j,k)+u2(3:i2+1,j,k)-u2(4:i2+2,j,k)) enddo enddo call mpi_sendrecv(u2(i2,1,1),1,m_vlon3d,m_e,1,scr1(i0,1,1),1,m_vlon3d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(u2(3,1,1),1,m_vlon3d,m_w,1,scr1(i1,1,1),1,m_vlon3d,m_e,1,m_cart,istat,ierr) do k=1,k1 do j=2,j1 scr(1,j,k)=scr(1,j,k)+iu(1,j,k)* (-scr1(i0,j,k)+u2(1,j,k)+u2(2,j,k)-u2(3,j,k)) scr(i1,j,k)=scr(i1,j,k)+iu(i1,j,k)*(-u2(i2,j,k)+u2(i1,j,k)+u2(i0,j,k)-scr1(i1,j,k)) !do i=1,i1 ! u(i,j,k)=o12*scr(i,j,k)*iu(i,j,k) !enddo u(1:i1,j,k)=o12*scr(1:i1,j,k)*iu(1:i1,j,k) enddo enddo do k=1,k1 !do j=js,jf ! do i=2,i1 ! scr(i,j,k)=6.*(v2(i,j,k)+v2(i,j+1,k)) ! enddo !enddo scr(2:i1,js:jf,k) = 6. * (v2(2:i1,js:jf,k)+v2(2:i1,js+1:jf+1,k)) !do j=2,j2 ! do i=2,i1 ! scr(i,j,k)=scr(i,j,k)-v2(i,j-1,k)+v2(i,j,k)+v2(i,j+1,k)-v2(i,j+2,k) ! enddo !enddo scr(2:i1,2:j2,k)=scr(2:i1,2:j2,k)-v2(2:i1,1:j2-1,k)+v2(2:i1,2:j2,k)+v2(2:i1,3:j2+1,k)-v2(2:i1,4:j2+2,k) enddo call mpi_sendrecv(v2(1,j2,1),1,m_vlat3d,m_n,1,scr1(1,j0,1),1,m_vlat3d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(v2(1,3,1),1,m_vlat3d,m_s,1,scr1(1,j1,1),1,m_vlat3d,m_n,1,m_cart,istat,ierr) !apply symmetry BC at North Pole !if (lopen==2) then ! call mpi_sendrecv(v2(1,j2,1),1,m_vlat3d,m_nd,1,scr1(1,j0,1),1,m_vlat3d,m_ns,1,m_cart,istat,ierr) ! call mpi_sendrecv(v2(1,3,1),1,m_vlat3d,m_ns,1,scr1(1,j1,1),1,m_vlat3d,m_nd,1,m_cart,istat,ierr) !endif do k=1,k1 ! skip loop for 2nd order scheme if (mylat .ne. 0) then do i=2,i1 scr(i,1,k)=scr(i,1,k)-scr1(i,j0,k)+v2(i,1,k)+v2(i,2,k)-v2(i,3,k) enddo endif !if (lopen==2) then ! do i=2,i1 ! scr(i,j1,k)=scr(i,j1,k)-v2(i,j2,k)+v2(i,j1,k)+v2(i,j0,k)-scr1(i,j1,k) ! enddo !else if (mylat .ne. ngy1) then do i=2,i1 scr(i,j1,k)=scr(i,j1,k)-v2(i,j2,k)+v2(i,j1,k)+v2(i,j0,k)-scr1(i,j1,k) enddo endif !endif !do j=js,jf ! do i=2,i1 ! v(i,j,k)=o12*scr(i,j,k)*iv(i,j,k) ! enddo !enddo v(2:i1,js:jf,k)=o12*scr(2:i1,js:jf,k)*iv(2:i1,js:jf,k) enddo ! the following non-porous lid approach may be replaced by a good estimate ! of e-p, as in medina and north pacifi! codes, as per dietrich, et al. ! e-p may be significant when there is significant lateral water vapor ! transport and fallout by translating hurricanes and typhoons. ! non-porous lid approach: ! maintain top layer salinity by exchange with layer 2, and ignore ! possible vortex stretching implied by non-zero e-p (w(i,j,1)). this ! simplifies global incompressibility treatment in domain decomposed mpi ! code. ! do 627 j=2,j1 ! do 627 i=2,i1 ! 627 w(i,j,1)=0. ! maintain surface climatological s in the long run ! 90 day nudging time scale toward surface climatological s !conserves total salt material assuming all wet points > one layer deep ! tmp1=1./(90.*daodt) ! tmp2=tmp1*odz(2)/odz(1) ! do 628 j=2,j1 ! do 628 i=2,i1 ! tmp=in(i,j,1)*(sclis(i,j,1,new)-s1(i,j,1)) ! s1(i,j,1)=s1(i,j,1)+tmp1*tmp ! 628 s1(i,j,2)=s1(i,j,2)-tmp2*tmp ! ==================================================== ! create non-divergent u,v,w for next step by ! "clearing out the divergence of the barotropic mode" ! ==================================================== ! use scr and scr to store changes for incompressibility do k=1,4 do j=1,j0 do i=1,i0 scr(i,j,k)=0. enddo enddo enddo ! --------------------------------------------------- ! iterate evp solver to get zero normal flow at shore ! and exactly non-divergent bt mode ! --------------------------------------------------- ! mxmask is maximum number of iterations allowed !check whether increasing from 8 to 9 eliminates noisy regions mxmask=9 ! note: it would be slightly more efficient to separate out the bt mode ! (as per dan wright's suggestion) if more than one iteration is ! performed, but this has not yet been implemented with river sources ! and precipitation call mpi_barrier(m_cart,ierr) timer(8)=mpi_wtime() itmask=0 652 itmask=itmask+1 ! eventually, we should use wright's bt model approach ! for efficiencfcy, but the bt mode's non-zero divergence ! requires special treatment in that case. ! zero out advection velocity over land tmp=0. do j=2,j1 do i=1,i1 u(i,j,1)=iu0(i,j)*u(i,j,1) enddo enddo do j=js,jf do i=2,i1 v(i,j,1)=iv0(i,j)*v(i,j,1) enddo enddo if (lfsrf .eq. 0) then do k=1,k1 temp1=1./odz(k) do j=2,j1 temp=ocs(j)*ody(j) do i=2,i1 w(i,j,k+1)=w(i,j,k)-((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*temp1 end do end do end do do j=2,j1 jj=j-1 !(j:cell, jj:size of s, easy to read) do i=2,i1 ii=i-1 !(i:cell, ii:size of s, easy to read) s(ii,jj)=w(i,j,kb(i,j)+1) end do end do elseif (lfsrf .eq. 1) then do k=k1,1,-1 temp1=1./odz(k) do j=2,j1 temp=ocs(j)*ody(j) do i=2,i1 w(i,j,k)=w(i,j,k+1)+((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*temp1 end do end do end do do j=2,j1 jj=j-1 !(j:cell, jj:size of s, easy to read) do i=2,i1 ii=i-1; !(i:cell, ii:size of s, easy to read) s(ii,jj)=-w(i,j,1) end do end do else if(myid.eq.0)write(*,*)'lfsrf is not set correctly' call mpi_finalize(ierr) stop endif ! ================================ ! the rebirth of bir ! something for m.s. to chew on!!! ! ================================ ! iter is the bir sweep number timer(9)=mpi_wtime() if (lsolver .eq. 1) then call rep(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie) ! call rep2(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie) elseif(lsolver .eq. 2) then #ifdef ALGO11 call p_bicgstab_le_algo11(ab,al,ac,ar,at,s,x,cb,cl,cc,cr,ct) #else call p_bicgstab_le(ab,al,ac,ar,at,s,x,cb,cl,cc,cr,ct,cgr,cgrh,cgp,cgv,cgs,cgt,cgph,cgsh) #endif endif ! call rep(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie) timer(9)=mpi_wtime()-timer(9) ! rigid-lid pressure adjustment !do j=1,j0 !do i=1,i0 ! p0(i,j)=p0(i,j)+odt*x(i,j) !enddo !enddo p0(1:i0, 1:j0) = p0(1:i0, 1:j0) + odt * x(1:i0, 1:j0) ! bt mode adjustment (from pressure gradient change) !do j=2,j1 !do i=1,i1 ! scr(i,j,1)=-(x(i+1,j)-x(i,j))*odx(j) !enddo !enddo do j=2,j1 scr(1:i1,j,1)=-(x(2:i1+1,j)-x(1:i1,j))*odx(j) enddo !do j=js,jf !do i=2,i1 ! scr(i,j,2)=-(x(i,j+1)-x(i,j))*odyv(j) !enddo !enddo do j=js,jf scr(2:i1,j,2)=-(x(2:i1,j+1)-x(2:i1,j))*odyv(j) enddo ! top layer !do j=2,j1 !do i=1,i1 ! scr(i,j,3)=scr(i,j,3)+scr(i,j,1) ! u(i,j,1)=iu0(i,j)*u(i,j,1)+scr(i,j,1) !enddo !enddo scr(1:i1, 2:j1, 3) = scr(1:i1, 2:j1, 3) + scr(1:i1, 2:j1, 1) u(1:i1, 2:j1, 1) = iu0(1:i1, 2:j1) * u(1:i1, 2:j1, 1) + scr(1:i1, 2:j1, 1) !do j=js,jf !do i=2,i1 ! scr(i,j,4)=scr(i,j,4)+scr(i,j,2) ! v(i,j,1)=iv0(i,j)*v(i,j,1)+scr(i,j,2) !enddo !enddo scr(2:i1, js:jf, 4) = scr(2:i1, js:jf, 4) + scr(2:i1, js:jf, 2) v(2:i1, js:jf, 1) = iv0(2:i1, js:jf) * v(2:i1, js:jf, 1) + scr(2:i1, js:jf, 2) ! remaining layers !do k=2,k1 ! do j=2,j1 ! do i=1,i1 ! u(i,j,k)=u(i,j,k)+scr(i,j,1)*iu(i,j,k) ! enddo ! enddo ! do j=js,jf ! do i=2,i1 ! v(i,j,k)=v(i,j,k)+scr(i,j,2)*iv(i,j,k) ! enddo ! enddo !enddo do k=2,k1 u(1:i1, 2:j1, k) = u(1:i1, 2:j1, k) + scr(1:i1, 2:j1, 1) * iu(1:i1, 2:j1, k) v(2:i1, js:jf, k) = v(2:i1, js:jf, k) + scr(2:i1, js:jf, 2) * iv(2:i1, js:jf, k) enddo !check convergence to zero flow over land tmp=0. do j=2,j1 do i=2,i2 tmp=max(tmp,(1-iu0(i,j))*abs(u(i,j,1))) enddo enddo do j=2,j2 do i=2,i1 tmp=max(tmp,(1-iv0(i,j))*abs(v(i,j,1))) enddo enddo tp=tmp call mpi_allreduce(tp,tmp,1,mpi_real8,mpi_max,m_cart,ierr) if(itf-it0.lt.itfday.or.mod(itf,itfday).eq.0) write(14,683) itf-it0,itmask,tmp,timer(9) !if (myid .eq. 0 .and. mod(itf,itfday).eq.0) write(*,683) itf-it0,itmask,tmp,timer(6) if (myid .eq. 0) write(*,683) itf-it0,itmask,tmp,timer(9) 683 format('@itf-it0===',i8,',itmask=',i3,', vmx on land=',f8.4,' cm/s, solver time =', f8.4,'sec') if (tmp.gt..0002 .and. itmask.lt.mxmask) go to 652 timer(8)=mpi_wtime()-timer(8) ! this completes the advanced time level advection velocity ! having exactly zero barotropic divergence and satisfying all ! kinematic and specified open boundary conditions if (lfsrf .eq. 0) then do k=1,k1 !(k:layer& top face for w, k+1:bottom face for w) tmp=1./odz(k) do j=2,j1 !(j:cell & north face for v) temp=ocs(j)*ody(j) do i=2,i1 !(i:cell & east face for u) w(i,j,k+1)=iw(i,j,k+1)*(w(i,j,k)-((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*tmp) end do end do end do elseif (lfsrf .eq. 1) then do j=2,j1 !(j:cell & north face for v) temp=ocs(j)*ody(j) !TS_NEW do i=2,i1 !(i:cell & east face for u) do k=kb(i,j),2,-1 tmp=1./odz(k) w(i,j,k)=(w(i,j,k+1)+((u(i,j,k)*iu(i,j,k)-u(i-1,j,k)*iu(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)*iv(i,j,k)-csv(j-1)*v(i,j-1,k)*iv(i,j-1,k))*temp)*tmp) end do k=1 tmp=1./odz(k) w(i,j,k)=(w(i,j,k+1)+((u(i,j,k)*iu0(i,j)-u(i-1,j,k)*iu0(i-1,j))*odx(j)+(csv(j)*v(i,j,k)*iv0(i,j)-csv(j-1)*v(i,j-1,k)*iv0(i,j-1))*temp)*tmp) end do end do endif do k=1,k1 do j=2,j1 do i=2,i1 u2(i,j,k)=12.*(u(i-1,j,k)+u(i,j,k)) v2(i,j,k)=12.*(v(i,j-1,k)+v(i,j,k)) enddo enddo enddo call mpi_sendrecv(u(i2,1,1),1,m_vlon3d,m_e,1,scr(i2,1,1),1,m_vlon3d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(u(2,1,1),1,m_vlon3d,m_w,1,scr(2,1,1),1,m_vlon3d,m_e,1,m_cart,istat,ierr) do k=1,k1 do j=2,j1 u2(2,j,k)=o24*in(2,j,k)*(u2(2,j,k)-scr(i2,j,k)+u(1,j,k)+u(2,j,k)-u(3,j,k)) u2(i1,j,k)=o24*in(i1,j,k)*(u2(i1,j,k)-u(i0-3,j,k)+u(i2,j,k)+u(i1,j,k)-scr(2,j,k)) enddo enddo do k=1,k1 do j=2,j1 do i=3,i2 u2(i,j,k)=o24*in(i,j,k)*(u2(i,j,k)-u(i-2,j,k)+u(i-1,j,k)+u(i,j,k)-u(i+1,j,k)) enddo enddo enddo call mpi_sendrecv(v(1,j2,1),1,m_vlat3d,m_n,1,scr(1,j2,1),1,m_vlat3d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(v(1,2,1),1,m_vlat3d,m_s,1,scr(1,2,1),1,m_vlat3d,m_n,1,m_cart,istat,ierr) !apply symmetry BC at North Pole !if (lopen==2) then ! call mpi_sendrecv(v(1,j2,1),1,m_vlat3d,m_nd,1,scr(1,2,1),1,m_vlat3d,m_ns,1,m_cart,istat,ierr) ! scr(:,2,:)=-scr(:,2,:) !endif do k=1,k1 do i=2,i1 do j=3,j2 v2(i,j,k)=v2(i,j,k)-v(i,j-2,k)+v(i,j-1,k)+v(i,j,k)-v(i,j+1,k) enddo if(mylat.ne.0) v2(i,2,k)=v2(i,2,k)-scr(i,j2,k)+v(i,1,k)+v(i,2,k)-v(i,3,k) !if(lopen==2) then ! v2(i,j1,k)=v2(i,j1,k)-v(i,j0-3,k)+v(i,j2,k)+v(i,j1,k)-scr(i,2,k) !else if(mylat.ne.ngy1) v2(i,j1,k)=v2(i,j1,k)-v(i,j0-3,k)+v(i,j2,k)+v(i,j1,k)-scr(i,2,k) !endif enddo enddo do k=1,k1 do j=js,j1 do i=2,i1 v2(i,j,k)=in(i,j,k)*o24*v2(i,j,k) enddo enddo enddo if (itf-it0.gt.itfday.and.mod(itf,itfday).ne.0) go to 712 !check incompressibility tmp=0. err=0. do k=1,k1 do j=2,j1 do i=2,i1 tmp1=(u(i,j,k)-u(i-1,j,k))*odx(j) tmp2=(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*ocs(j)*ody(j) tmp3=(w(i,j,k+1)-w(i,j,k))*odz(k) tmp=tmp+max(abs(tmp1),abs(tmp2),abs(tmp3))*in(i,j,k) err=err+abs(tmp1+tmp2+tmp3)*in(i,j,k) enddo enddo enddo tep=err call mpi_allreduce(tep,err,1,mpi_real8,mpi_sum,m_cart,ierr) tep=tmp call mpi_allreduce(tep,tmp,1,mpi_real8,mpi_sum,m_cart,ierr) err=err/tmp write(14,711) err call flush(14) 711 format(' *** normalized mean incompressibility error = ',1pe9.2) 712 psm=0. do j=2,j1 do i=2,i1 psm=psm+p0(i,j) enddo enddo tmp=psm call mpi_allreduce(tmp,psm,1,mpi_real8,mpi_sum,m_cart,ierr) psm=psm/dble(i2t*j2t) do j=2,j1 do i=2,i1 p0(i,j)=p0(i,j)-psm enddo enddo call mpi_barrier(m_cart,ierr) timer(10)=mpi_wtime() call mpi_exch_3d_r8(u2) if (lopen==2) then call mpi_exch_3d_r8_vsymtry(v2) else call mpi_exch_3d_r8(v2) endif call mpi_exch_3d_r8(t2) call mpi_exch_3d_r8(s2) call mpi_exch_3d_r8(c2) timer(10)=mpi_wtime()-timer(10) if (.false.) then tt0=it0/daodt n=tt0/360.+.0001 tmp=tt0/360. tt2=(tmp-n)*360.+.1 tt1=tt0/360.+.0001+1 write(fn1,'(i3.3,a1,i3.3)'),tt1,'_',tt2 do i = 1,i0,2 do j = 1,j0,2 if (ydeg(j).ge.30.1.and.ydeg(j).le.40.1 .and.xdeg(i).ge.(360.-80.).and.xdeg(i).le.(360.-72.)) then open(1005,file='ts/field1_'//fn1//'_'//grd,form='unformatted') do k=1,10 write(1005) ulf(i,j,k),vlf(i,j,k),p(i,j,k),tlf(i,j,k) enddo if(itf.eq.it0) then print*,'field i,j,myid',myid,grd,i,j,ydeg(j),xdeg(i),ulf(i,j,4) endif endif if (ydeg(j).ge.36.1.and.ydeg(j).le.42.1 .and.xdeg(i).ge.(360.-70.).and.xdeg(i).le.(360.-60.)) then open(1006,file='ts/field2_'//fn1//'_'//grd,form='unformatted') do k=1,10 write(1006) ulf(i,j,k),vlf(i,j,k),p(i,j,k),tlf(i,j,k) enddo if(itf.eq.it0) then print*,'field i,j,myid',myid,grd,i,j,ydeg(j),xdeg(i),ulf(i,j,4) endif endif enddo enddo endif ! totalt=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt ! call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT7',TOTALT/TOTALVOL !20210904 ! ============================= ! heat flux from ice-formation ! ============================= if (liceform) then timer(13)=mpi_wtime() call ice_formation timer(13)=mpi_wtime()-timer(13) !if (lactive_ice) call ice_flx_to_coupler !WRITE(*,*)'[fsglo]check_time_flag(ice_cpl_flag)=',check_time_flag(ice_cpl_flag),nsteps_run ! if ( check_time_flag(ice_cpl_flag) ) then ! call ice_flx_to_coupler ! endif end if ! rawf=0.53 ! fltw=0.1 ! modified filter scheme call modified_filter !totalt=sum(t1(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) !dtmp=totalt !call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) !write(40,*) 't1_2',TOTALT/TOTALVOL !TS if(.false.)then ! ---------------------------------------------- ! relax swamp layer and masked water to interior ! to reduce seepage effect ! ---------------------------------------------- call mpi_barrier(m_cart,ierr) call mpi_exch_2d_r8(t1(1,1,1)) call mpi_exch_2d_r8(s1(1,1,1)) call mpi_exch_2d_r8(c1(1,1,1)) do j=2,j1 do i=2,i1 tmp=.25*(1-in(i,j,1)) scr(i,j,1)=in(i,j,1)*t1(i,j,1)+tmp*(t1(i-1,j,1)+t1(i+1,j,1)+t1(i,j+1,1)+t1(i,j-1,1)) scr(i,j,2)=in(i,j,1)*s1(i,j,1)+tmp*(s1(i-1,j,1)+s1(i+1,j,1)+s1(i,j+1,1)+s1(i,j-1,1)) scr(i,j,13)=in(i,j,1)*c1(i,j,1)+tmp*(c1(i-1,j,1)+c1(i+1,j,1)+c1(i,j+1,1)+c1(i,j-1,1)) enddo enddo !do j=2,j1 !do i=2,i1 ! t1(i,j,1)=scr(i,j,1) ! s1(i,j,1)=scr(i,j,2) ! c1(i,j,1)=scr(i,j,13) !enddo !enddo t1(2:i1, 2:j1, 1) = scr(2:i1, 2:j1, 1) s1(2:i1, 2:j1, 1) = scr(2:i1, 2:j1, 2) c1(2:i1, 2:j1, 1) = scr(2:i1, 2:j1, 13) endif ! totalt=sum(t2(2:i1,2:j1,:)*volume(2:i1,2:j1,:)) ! dtmp=totalt ! call mpi_allreduce(dtmp,totalt,1,mpi_real8,mpi_sum,m_cart,ierr) ! write(23,*) 'meanT8',TOTALT/TOTALVOL end subroutine subroutine set_windmix ! open(8,file='prepglo/data/windmix'//grd,form='unformatted') ! read(8) vbk,hbk ! close(8) call rdnc_windmix if (lkppmix .and. ltidal_mixing) then call rdnc_tidal !----------------------------------------------------------------------- ! convert TIDAL_ENERGY_FLUX from W/m^2 to gr/s^3 and apply mask !----------------------------------------------------------------------- tidal_energy_flux(:,:)=1000.d0*in(:,:,1)*tidal_energy_flux(:,:) endif end subroutine subroutine set_rundata use time_management, only: iyear,imonth ! open(9,file='prepglo/data/rundata'//grd,form='unformatted') ! read(9) case,dscrib,dt,tlz,b,g,dm0,de0,dmz0,f,tanphi,rzmx, & ! drag,fltw,daodt,odt,prn,orzmx,ofltw,y,yv,yvdeg,ydeg,cs,csv, & ! ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv,ktrm,lrstrt,mxit,isav, & ! mxsav,lopen,lmovi,mvi,lfsrf,lsolver,xdeg,lprecon,rawf ! close(9) call rdnc_rundata if (myid .eq. 0) then ! open(79,file='namelist.txt',form='formatted') ! read(79,*)lrstrt ! read(79,*)daodt ! read(79,*)irun ! close(79) ! mxit=irun*daodt*365 dt=2.d0*86400.d0/daodt endif call mpi_bcast(lrstrt,1,mpi_int,0,m_cart,ierr) call mpi_bcast(daodt,1,mpi_real8,0,m_cart,ierr) call mpi_bcast(mxit,1,mpi_int,0,m_cart,ierr) call mpi_bcast(dt,1,mpi_int,0,m_cart,ierr) itfday=daodt+.0001 !initialize the month and the relevent time avg. variable nmon=ini_clim if (lrstrt==0.and.iyear==1.and.imonth==1) then nopt=dmonth(nmon)-1 else nopt=dmonth(nmon) endif wopt=1./dble(nopt*itfday) do nm=1,12 ttmon(nm)=sum(dmonth(1:nm)) epday(nm)= ttmon(nm) - 1.d0/daodt if(myid==0) write(*,*)'epday=',epday(nm),nm enddo !modify the meridional boundary points index for lopen !lopen==2 for symmetry BC at North Pole if (lopen==2 .and. mylat==ngy1) jf=j0-1 !conversion of qavg (deg c per model time step) to watts/m-m ! use factors for deg ! per year and g/kg/yr watts=360.*24.*3600.*0.00134/(dt*odz(1)) ! global salinity budget weighting factor end subroutine subroutine set_zkb ! -------------------------------------------------------------------- ! metrics, z-levels(z) and bathymetry (logical depth array (kb)) data ! -------------------------------------------------------------------- ! open(60,file='prepglo/data/zkb'//grd,form='unformatted') ! read(60) tlz,z,odz,odzw,kb,in,iu,iv,iw,iu0,iv0 call rdnc_zkb end subroutine subroutine set_surfwinds ! open(72,file='prepglo/data/winds'//grd,form='unformatted') call rdnc_winds do n=1,12 ! read wind stress (dynes/cm-cm), temperature (deg c), salinity(ppt), ! q, wavg(mm/month) ! read(72) ((taux(i,j,n),i=1,i2),j=1,j2),((tauy(i,j,n),i=1,i2),j=1,j2) ! screen out land sources and convert from newtons/m-m to dynes/cm-cm ! intensifies the mixing le nocn=0 do j=1,j2 do i=1,i2 nocn=nocn+in(i+1,j+1,1) taux(i,j,n)=10.*in(i+1,j+1,1)*taux(i,j,n) tauy(i,j,n)=10.*in(i+1,j+1,1)*tauy(i,j,n) enddo enddo enddo nt=nocn call mpi_allreduce(nt,nocn,1,mpi_integer,mpi_sum,m_cart,ierr) end subroutine subroutine set_solarad call rdnc_rad do n=1,12 do j=1,j2 do i=1,i2 qsw(i,j,n)=in(i+1,j+1,1)*qsw(i,j,n) enddo enddo enddo end subroutine subroutine set_area_volume real*8 dtmp do k=1,k1 a(k)=0.d0 do j=2,j1 dtmp=dx(j)*dy(j) do i=2,i1 a(k)=a(k)+in(i,j,k)*dtmp enddo enddo dtmp=a(k) call mpi_allreduce(dtmp,a(k),1,mpi_real8,mpi_sum,m_cart,ierr) enddo volume=0.d0 do k=1,k1 do j=2,j1 dtmp=dx(j)*dy(j)/odz(k) do i=2,i1 volume(i,j,k)=dtmp*in(i,j,k) enddo enddo enddo dtmp=sum(volume(2:i1,2:j1,:)) call mpi_allreduce(dtmp,totalvol,1,mpi_real8,mpi_sum,m_cart,ierr) dtmp=sum(volume(2:i1,2:j1,1:5)) call mpi_allreduce(dtmp,totalvol_sfc,1,mpi_real8,mpi_sum,m_cart,ierr) dtmp=sum(volume(2:i1,2:j1,1)) call mpi_allreduce(dtmp,totalvol1,1,mpi_real8,mpi_sum,m_cart,ierr) end subroutine subroutine set_ts_nudging ! simple equatorial and high latitude sponge layers are needed because ! storm-generated penetrative convection in high lats and inappropriate ! hydrostati! equatorial dynamics that cannot emulate effe! of stacked ! equatorial vortices on t,s near equator; the effe! on momentum are ! addressed by new eddy viscosity formulation based on the time ! variability of vertical velocity which is big: near the equator (due to ! small coriolis parameter); in high lats and southern ocean due to small ! stratification there; and along fronts everywhere due to instabilities ! ====================================================================== ! jacc,jeq,jar! refer to cell centers. assume southernmost and ! northernmost gridlines are at the same lats as used at 1/4 deg lon ! resolution, which is achieved by proper inputs, defined in ! prepglo/data/metgen. ! 1. keep the same southernmost gridline lat (y0deg=yv(1)=-69.9377 deg) ! as used with 1/4 deg resolution (dxmnut=15.) ! 2. set the new j2 value equal to the old one (792) multiplied by the ! ratio of the old "dxdeg" value to the new one, ! this puts a gridline exactly on the equator and puts the northernmost ! gridline, yv(j1), at 69.7654 deg lat; 1/4 deg res (dxmnut=15.), gives ! yv(398)= 0.0 deg (on the equator). ! ====================================================================== ! edit: for now, T&S nudging intensity are constant across lats. ! --------------------- ! setting nudge matrix ! --------------------- real*8 n_tmp t_nudstr=max(1./(1080.d0*daodt),0.) s_nudstr=max(1./(1080.d0*daodt),0.) scr=0. yscale=10. rate=(792./j2t)/yscale rate3=3.*rate ! s n_tmp=-1. tmpp_s=max(1./(n_tmp*daodt),0.) ! t n_tmp=-1. tmpp_t=max(1./(n_tmp*daodt),0.) ! s,t coast n_tmp=-1. tmpp_sc=max(1./(n_tmp*daodt),0.) n_tmp=-1. tmpp_tc=max(1./(n_tmp*daodt),0.) ! s,t ac! arc n_tmp=-1. !5 tmpp_shs=max(1./(n_tmp*daodt),0.) n_tmp=-1. !5 , 0926 meeting tmpp_ths=max(1./(n_tmp*daodt),0.) ! hmean ac! arc ! n_tmp=365. n_tmp=-1. tmpp_hi=max(1./(n_tmp*daodt),0.) ! hmean eq n_tmp=-1. tmpp_ei=max(1./(n_tmp*daodt),0.) ! hmean other n_tmp=-1. tmpp_o=max(1./(n_tmp*daodt),0.) ! hmean coast n_tmp=-1. tmpp_oc=max(1./(n_tmp*daodt),0.) acc_deg=-65. acc_rag=acc_deg+2. arc_deg=70. arc_rag=arc_deg-2. do j=2,j1 do i=2,i1 ! s & t if(ydeg(j).ge.arc_deg.or.ydeg(j).le.acc_deg) then t_nudge(i,j)=tmpp_ths s_nudge(i,j)=tmpp_shs elseif((ydeg(j).ge.arc_rag)) then t_nudge(i,j)=(tmpp_ths-tmpp_t)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_t s_nudge(i,j)=(tmpp_shs-tmpp_s)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_s elseif((ydeg(j).le.acc_rag)) then t_nudge(i,j)=(tmpp_ths-tmpp_t)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_t s_nudge(i,j)=(tmpp_shs-tmpp_s)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_s else t_nudge(i,j)=tmpp_t s_nudge(i,j)=tmpp_s endif if(kb(i,j).le.1.or.kb(i-1,j).le.1.or.kb(i+1,j).le.1.or.kb(i,j-1).le.1.or.kb(i,j+1).le.1) then t_nudge(i,j)=tmpp_tc s_nudge(i,j)=tmpp_sc endif ! hmean do k=1,k0 if(ydeg(j).ge.arc_deg.or.ydeg(j).le.acc_deg) then hmean_nudge(i,j,k)=tmpp_hi elseif((ydeg(j).ge.arc_rag)) then hmean_nudge(i,j,k)=(tmpp_hi-tmpp_o)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_o elseif((ydeg(j).le.acc_rag)) then hmean_nudge(i,j,k)=(tmpp_hi-tmpp_o)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_o elseif(ydeg(j).ge.-5..and.ydeg(j).le.5.) then hmean_nudge(i,j,k)=(tmpp_ei-tmpp_o)*exp(-rate3*abs((ydeg(j)-0.)/(ydeg(j)-ydeg(j-1))))+tmpp_o else hmean_nudge(i,j,k)=tmpp_o endif if(kb(i,j).le.1.or.kb(i-1,j).le.1.or.kb(i+1,j).le.1.or.kb(i,j-1).le.1.or.kb(i,j+1).le.1) then hmean_nudge(i,j,k)=tmpp_oc endif enddo enddo enddo end subroutine subroutine set_lateral_mixing ! increase lateral mixing near equator and in high latitudes ! big northern latitude dissipation fixes initialization data problem it=mylon*i2 jt=mylat*j2 dxx=1.0 dyy=1.0 do k=1,k1 do j=1,j0 do i=1,i0 temp=1.e7*(exp(-.001*(2-j-jt)**2)+exp(-.001*(j1t-jt-j)**2)) ! north & south poles scr(i,j,k)=2.e7*exp(-.001*((ydeg(j)-0)/dyy)**2)+temp !+.5e6*exp(-.01*(((xdeg(i)-109.5))**2+((ydeg(j)-15.5))**2)) &! warm bias !+.5e6*exp(-.01*(((xdeg(i)-114.5))**2+((ydeg(j)+ 9.5))**2)) &! warm bias !+2.e6*exp(-.01*(((xdeg(i)-300.9))**2+((ydeg(j)-45.11))**2))&! cold bias (west nova scotia) !+2.e6*exp(-.01*(((xdeg(i)-298.7))**2+((ydeg(j)-43.68))**2))&! cold bias (sw nova scotia) !+2.e6*exp(-.01*(((xdeg(i)-276.2))**2+((ydeg(j)-56.67))**2))&! cold bias (south hudson bay) !+2.e6*exp(-.01*(((xdeg(i)-181.7))**2+((ydeg(j)-64.22))**2))&! cold bias (gulf of anadyr) !+2.e6*exp(-.01*(((xdeg(i)-164.8))**2+((ydeg(j)-58.85))**2))&! cold bias (ne kamchatka) !+1.e6*exp(-.01*(((xdeg(i)-142.3))**2+((ydeg(j)-56.67))**2))&! cold bias (sakhalin gulf) !+1.e6*exp(-.01*(((xdeg(i)-143.4))**2+((ydeg(j)-45.11))**2))&! cold bias (aniva bay) !+1.e6*exp(-.01*(((xdeg(i)-135.6))**2+((ydeg(j)-42.21))**2)) ! cold bias (east japan sea) enddo enddo enddo ! note: the following may need some adjustment in the mpi code ! increase lateral eddy viscosity near wet/dry point interfaces do k=1,k1 do i=1,i0 do j=1,j0 dmx(i,j,k)=dm0+scr(i,j,k) dmy(i,j,k)=dmx(i,j,k) enddo enddo enddo end subroutine subroutine set_extra_vertical_eddy_viscosity ! xtra vertical eddy viscosity do jj=2,j1 j=min(jj,j2) j=max(j,3) do ii=2,i1 i=min(ii,i2) i=max(i,3) ! bottom slope dzdx=(z(2*kb(i+1,j))-z(2*kb(i-1,j)))*odx(j) dzdy=(z(2*kb(i,j+1))-z(2*kb(i,j-1)))*ody(j) zbot=z(2*kb(i,j)+1) slope=sqrt(dzdx**2+dzdy**2) do k=1,k2 ! regionally increased vertical eddy viscosity ! locally increase add array to get 90-day time scale for 1 layer scale ! tscale=90.*24.*3600. ! after day 100 yr 3 confine xtra equatorial (j=398) mixing to top 1000m tmp=min(10.,2.e-4*z(2*kb(i,j)+1)) ! evadd=.5*exp(-tmp-.01*(jt+j-j2t/2-2)**2) evadd=.5*exp(-tmp-.01*((ydeg(j)-0)/(ydeg(j+1)-ydeg(j)))**2) ! eddy viscosity is big above steeply sloping bottom locations ! near-bottom factor ! temp=2.e-9*(zbot-z(2*k+1))**2 ! decrease vertical mixing scale at start of yr 3 to allow free ! passage into caribbean temp=1.e-8*(zbot-z(2*k+1))**2 ! vertical mix time scale=o(dz*2/ev); for time scale=1 day and dz=10m, ! ev=o(10**6/10**5)=10. tmp=min(5.,5.*exp(-temp)*slope) ! biggest eddy viscosity in deep water having steeply sloping bottom ! where breaking big amplitude internal waves are most likely to occur ! due to low subthermocline stratification tmp=tmp+vbk(k)+evadd ! explicit stability limit and ev normalization done in pp82glo add(ii-1,jj-1,k)=tmp add=0.d0 enddo enddo enddo end subroutine subroutine set_seasonal_data use time_management, only: iyear,imonth real*8 tmpp ! open(70,file='prepglo/data/annualevitus'//grd,form='unformatted') ! rewind 70 ! seasonal boundary conditions do nm=1,12 if(myid.eq.0)write(*,15) nm 15 format('month',i3,' in bounds') ! read in monthly t,s boundary climatologies ! read(70) sclim,tclim call rdnc_levitus(nm) do k=1,k1 do j=1,j0 do i=1,i0 s2(i,j,k)=s2(i,j,k)+sclim(i,j,k)*dmonth(nm)/365.d0 t2(i,j,k)=t2(i,j,k)+tclim(i,j,k)*dmonth(nm)/365.d0 !xudog,1109 !if(k.eq.1) sclim(i,j,k)=max(33.,sclim(i,j,k)) if(ini_clim.eq.13)then s1(i,j,k)=s1(i,j,k)+sclim(i,j,k)*dmonth(nm)/365.d0 t1(i,j,k)=t1(i,j,k)+tclim(i,j,k)*dmonth(nm)/365.d0 else if(ini_clim.eq.nm)then s1(i,j,k)=sclim(i,j,k) t1(i,j,k)=tclim(i,j,k) endif endif enddo enddo enddo do j=1,j0 do i=1,i0 do k=1,k1 sclis(i,j,k,nm)=sclim(i,j,k) tclis(i,j,k,nm)=tclim(i,j,k) enddo enddo enddo enddo ! close(70) do k=1,k1 do j=1,j0 do i=1,i0 sclim(i,j,k)=s2(i,j,k) tclim(i,j,k)=t2(i,j,k) s2(i,j,k)=0. t2(i,j,k)=0. enddo enddo enddo if (lrstrt.eq.0)then if(myid.eq.0)write(*,*)'write ncf initial' call write_init2nc(sclim,tclim) do k=1,k1 do j=1,j0 do i=1,i0 slf(i,j,k)=s1(i,j,k) tlf(i,j,k)=t1(i,j,k) s2(i,j,k)=s1(i,j,k) t2(i,j,k)=t1(i,j,k) !yc c1(i,j,k)=0. !t1(i,j,k) !assume c1 is same as temp clf(i,j,k)=c1(i,j,k) c2(i,j,k)=c1(i,j,k) !yc enddo enddo enddo write(14,1)case,dscrib,lrstrt,mxit,i0,j0,k0,ktrm,dt,fltw,tlz,g,dm0,de0,dmz0,rzmx,drag 1 format('case ',a5,' (',a63,')'/'----------'/ & ' lrstrt mxit i0 j0 k0 ktrm'/2i9,4i5// & 8x,'dt',6x,'fltw',7x,'tlz'/1p,3(1x,e9.2)// & 7x,'g',7x,'dm0',7x,'de0',6x,'dmz0', & 6x,'rzmx',6x,'drag'/1p,7(1x,e9.2)) ! initialize qavg and wavg !n=30*daodt temp=z(3)/(365.*24.*3600.) do m=1,12 nsombo(m)=1 !EPmodified,20171031 mp=m+1 if (m.eq.12) mp=1 do j=1,j2 do i=1,i2 ! this initialization is arbitrary and optional qavg(i,j,m)=in(i+1,j+1,1)*(tclis(i+1,j+1,1,mp)-tclis(i+1,j+1,1,m))*t_nudstr wavg(i,j,m)=in(i+1,j+1,1)*(sclis(i+1,j+1,1,mp)-sclis(i+1,j+1,1,m))*s_nudstr enddo enddo enddo !Initialize qavg and wavg temp array. EPmodified,20171031 qtmp=qavg wtmp=wavg ! depth ! read(60) ((scr(i,j,1),i=1,i0),j=1,j0) ! close(60) call rdnc_zkb2 do j=2,j1 do i=2,i1 ! convert to kilometers scr(i,j,1)=1.e-5*scr(i,j,1) enddo enddo else ! ----------------- ! read restart data ! ----------------- ! av is for long term averages. avg is for nabcs 5-day averages. if (myid.eq.0) print*,'read rerun file',myid,itf call read_rerun !derive nmonth for historical time-avg. output and scheme time management. syng_days=mod(days,365.d0) !find current month dtmp=syng_days nld=0 if(dtmp.eq.0) nld=12 do while(dtmp.gt.0) nld=nld+1 dtmp=dtmp-dmonth(nld) enddo if(nld.eq.0) then dtmp=-1.d0 nld=1 new=2 elseif (nld.eq.12) then new=1 else new=nld+1 endif fold=-1.d0*dtmp/dmonth(nld) fnew=1.d0-fold nmon=nld if (myid==0) write(*,*) 'ratio',syng_days,fnew,new,fold,nld !hist time avg. variables if (lrstrt==0.and.iyear==1.and.imonth==1) then nopt=dmonth(nmon)-1 else nopt=dmonth(nmon) endif wopt=1./dble(nopt*itfday) call mpi_exch_3d_r8(u) call mpi_exch_2d_r8(x) if (lopen==2) then call mpi_exch_3d_r8_vsymtry(v) else call mpi_exch_3d_r8(v) endif if(myid.eq.0)write(*,*)'read ncf initial' call read_initnc(sclim,tclim) call mpi_exch_3d_r8(tclim) call mpi_exch_3d_r8(sclim) n=days*daodt+.5 daodt_new=daodt daodt_old=itf/days if (n/=itf) then ! tmp=(dtnew/dtold) tmp=daodt_old/daodt_new if (myid .eq. 0)then write(*,57) tmp write(14,57) tmp 57 format('dtnew/dtold=',f6.3) endif do nm=1,12 do j=1,j2 do i=1,i2 qavg(i,j,nm)=tmp*qavg(i,j,nm) wavg(i,j,nm)=tmp*wavg(i,j,nm) qtmp(i,j,nm)=tmp*qtmp(i,j,nm) wtmp(i,j,nm)=tmp*wtmp(i,j,nm) enddo enddo enddo if (myid .eq. 0)then write(*,59) itf,n write(14,59) itf,n 59 format(/'time step # changed from',i7,' to',i7,' at restart to reflect dt change') endif itf=n endif if(myid.eq.0)write(*,61) itf,days,daodt_new,daodt_old 61 format('itf,days,daodt_new,daodt_old',i8,3f8.1) n=days+.1 if (mod(int(days+.01),1800).eq.0.) av=0. itf=days*daodt+0.5 ! 82 if (days.gt.460..and.mod(itf,360*itfday).ne.0) go to 85 ! wavg & qavg include advection effe! ! wavg is surface freshwater sources/sinks (e-p) consistent with model ! and climatological data ! nsombo(m)=number of ensemble-averaging years ! assumed already included in qavg and wa[vg ! however, we weight each new year equally to the previously completed ! year's ensemble average, until year 11, thus given more weight to ! recent results and thus acclerating convergence to a multi-year ensemble ! average. if (n<=10*365) then ! keep old avg but speed convergence to new avg do m=1,12 nsombo(m)=1 !EPmodified enddo endif if(myid .eq. 0)then write(*,90) itf,days,nsombo write(14,90) itf,days,nsombo 90 format(/'restart data read for time step ',i9,', day ',f8.2/'nsombo:',12i3) endif endif if (mylat .eq. 0) then do k=1,k1 do i=1,i0 t1(i,1,k)=tclim(i,1,k) t2(i,1,k)=tclim(i,1,k) tlf(i,1,k)=tclim(i,1,k) s1(i,1,k)=sclim(i,1,k) s2(i,1,k)=sclim(i,1,k) slf(i,1,k)=sclim(i,1,k) enddo enddo endif if (mylat .eq. ngy1) then do k=1,k1 do i=1,i0 t1(i,j0,k)=tclim(i,j0,k) t2(i,j0,k)=tclim(i,j0,k) tlf(i,j0,k)=tclim(i,j0,k) s1(i,j0,k)=sclim(i,j0,k) s2(i,j0,k)=sclim(i,j0,k) slf(i,j0,k)=sclim(i,j0,k) enddo enddo endif !call rdnc_STF end subroutine subroutine eliminate_spurious_values if(lrstrt.eq.0)then ! eliminate spurious values over land ! nit=2001 nit=1 do n=1,nit do k=1,2 tmin=1.e8 tmax=-1.e8 smin=1.e8 smax=-1.e8 do j=2,j1 do i=2,i1 if (in(i,j,k)==0) then scr(i,j,1)=.25*(s1(i+1,j,k)+s1(i-1,j,k)+s1(i,j-1,k)+s1(i,j+1,k)) smin=min(smin,scr(i,j,1)) smax=max(smax,scr(i,j,1)) scr(i,j,2)=.25*(t1(i+1,j,k)+t1(i-1,j,k)+t1(i,j-1,k)+t1(i,j+1,k)) tmin=min(tmin,scr(i,j,2)) tmax=max(tmax,scr(i,j,2)) endif enddo enddo do j=2,j1 do i=2,i1 if (in(i,j,k)==0) then s1(i,j,k)=scr(i,j,1) t1(i,j,k)=scr(i,j,2) endif enddo enddo call mpi_exch_2d_r8(t1(1,1,k)) call mpi_exch_2d_r8(s1(1,1,k)) if(myid.eq.0)then if(mod(n,100).eq.1)write(*,2) n,k,tmin,tmax,smin,smax 2 format('n,k,tmin,tmax,smin,smax:',i5,i3,4f7.3) endif enddo enddo endif end subroutine subroutine fw_conserv !global salt conservation epcomp= 0.d0 evap = 0.d0 dtmp = dt*odz(1)*1000.d0 !TSb dtmp = dtmp/2.0d0 do j=2,j1; jj=j-1; temp=dx(j)*dy(j) do I=2,i1; ii=i-1; epsrc(i,j)=temp*(wavg(ii,jj,nld)+snud_tmp(ii,jj)) + & dtmp*temp*STF(ii+1,jj+1,2) evap=evap+in(i,j,1)*epsrc(i,j) end do end do tmp=evap call mpi_allreduce(tmp,evap,1,mpi_real8,mpi_sum,m_cart,ierr) epcomp=evap/a(1) end subroutine subroutine ts_nudging_for_watermass ! To emulate watermass transformation effects on THC, nudge toward ! eqatorial as well as hi and lo lat climatology. We do not directly ! address watermass tranformations closer to poles, because their most ! important effects are seen in the central ACC and Arctic Circle ! regions, but important wind forced dynamics are included closer to ! poles. This is valid for simulating present climate and model ! validation, but NOT for simulating climate change, which equires ! ocean coupling to full atmospheric and ice dynamics, and much more ! complex physics. ! NOTE: This gives internal effective heat and salt material "sources" ! not in the real ocean. However, this should give model climatology ! close to the real ocean with the "sources" magnitudes being a metric ! of model quality; they should decrease with increased resolution and ! better model physics and subgridscale mixing emulation. Thus, the time ! averaged source magnitudes are a good indicator of model applicability ! to climate research. Indeed, the following may eventually include time ! averaged "sources" to decrease damping of transients (analagous to the ! non-damping surface sources of Dietrich, et al., 2004), but are not ! simply surface boundary conditions needed for ocean-only modeling (may ! be derived by coupling "perfect" ocean and atmospheri! models). They ! also noted that the "sources" may be used to improve vertical mixing ! parameterizations needed to keep imperfect models on track with ! observed ocean climate. ! ====================================================================== ! Later, check sensitivity to expanding nudging to bands of polar and ! eqatorial lats to invoke stronger watermass formation. ! Loop 95 replaces hardwired multi-region approach used by original ! coding before calling FSGLO. evap_hmean=0.d0 do k=2,k1 do j=2,j1 dv=dx(j)*dy(j)/odz(k) do i=2,i1 if(in(i,j,k)/=0) then t2(i,j,k)=(1.-hmean_nudge(i,j,k))*t2(i,j,k)+hmean_nudge(i,j,k)*tclim(i,j,k) s2(i,j,k)=(1.-hmean_nudge(i,j,k))*s2(i,j,k)+hmean_nudge(i,j,k)*sclim(i,j,k) !residual salinity that have to be subtracted later !to maintain avg.salinity ( evap_hmean in (cm^3/s) ) !ds=hmean_nudge(i,j,k)*(sclim(i,j,k)-s2(i,j,k)) !evap_hmean=evap_hmean + dv*ds/((s2(i,j,k)-ds)*dt) endif enddo enddo enddo do j=2,j1 dv=dx(j)*dy(j)/odz(1) do i=2,i1 t2(i,j,1)=t2(i,j,1)+t_nudge(i,j)*in(i,j,1)*in(i,j,2)*& (fold*tclis(i,j,1,nld)+fnew*tclis(i,j,1,new)-t2(i,j,1)) !20210329 move to set_STF ! s2(i,j,1)=s2(i,j,1)+s_nudge(i,j)*in(i,j,1)*in(i,j,2)*& ! (fold*sclis(i,j,1,nld)+fnew*sclis(i,j,1,new)-s2(i,j,1)) !residual salinity that have to be subtracted later !to maintain avg.salinity !ds=s_nudge(i,j)*in(i,j,1)*in(i,j,2)*(fold*sclis(i,j,1,nld)+fnew*sclis(i,j,1,new)-s2(i,j,1)) !evap_hmean=evap_hmean + dv*ds/((s2(i,j,1)-ds)*dt) enddo enddo end subroutine subroutine solve_hydrostatic_equation real*8 tmpd,sltd,pisd,p01,p02 dimension rhomin(k1) real*8, dimension(i2,j2,k1) :: wface ! ====================================================== ! Hydrostatic equation with static balance subracted out ! ====================================================== if (.not.lmkcof_eos) then do k=1,k1 temp=1.e8 do j=2,j1 do i=2,i1 tmpd=max(-1.8,t2(i,j,k)) sltd=min(41.,s2(i,j,k)) sltd=max(0.d0,sltd) pisd=.001*z(2*k) ! ======================= ! dan wright full e.o.s. ! rho(i,j,k)=r(tmpd,sltd,pisd) ! function r(t,s,p) ! ======================= pp02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd) pp01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd rho(i,j,k)=pp01/(pp02+0.7028423*pp01) ! determine global min of rho in all layers if(in(i,j,k)/=0) then temp=min(dble(temp),rho(i,j,k)) endif enddo enddo tmp=temp call mpi_barrier(m_cart,ierr) call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_min,m_cart,ierr) rhomin(k)=temp enddo else if ((lmkcof_eos) .and. (.not. eqstat_eos)) then !eos do 97 k=1,k1 dtemp=1.d8 pisd=1.d-3*z(2*k) do 96 j=1,j0 do 96 i=1,i0 if(in(i,j,k).eq.0) goto 96 tmpd=MAX(-1.8d0,t2(i,j,k)) tmpd=MIN(40.d0,tmpd) sltd=MIN(41.d0,s2(i,j,k)) sltd=MAX(0.d0,sltd) call make_local_coefficients & (tmpd,sltd,pisd,ccrho,ccTP,ccS,ccTPTP,ccSTP) rho(I,J,K)=local_density(tmpd,sltd,ccrho,ccTP,ccS,ccTPTP,ccSTP) & *1.d-3+1.d0 if(in(i,j,k).eq.0) goto 96 dtemp=MIN(dtemp,rho(i,j,k)) 96 CONTINUE dtmp=dtemp CALL MPI_BARRIER(M_CART,IERR) TIMER(9)=MPI_WTIME() CALL MPI_ALLREDUCE(DTMP,DTEMP,1,MPI_REAL8,MPI_MIN,M_CART,IERR) TIMER(3)=TIMER(3)+MPI_WTIME()-TIMER(9) rhomin(k)=dtemp 97 CONTINUE else ! eqstat_eos do k=1,k1 dtemp=1.d8 call state(k, k, t2(:,:,k), s2(:,:,k)*ppt_to_salt, rhofull=rho(:,:,k)) do j=1,j0 do i=1,i0 if(in(i,j,k)/=0) then dtemp=MIN(dtemp,rho(i,j,k)) endif enddo enddo dtmp=dtemp CALL MPI_BARRIER(M_CART,IERR) TIMER(9)=MPI_WTIME() CALL MPI_ALLREDUCE(DTMP,DTEMP,1,MPI_REAL8,MPI_MIN,M_CART,IERR) TIMER(3)=TIMER(3)+MPI_WTIME()-TIMER(9) rhomin(k)=dtemp enddo end if !eos ! Subtract off global min of rho (in all layers, but it is always ! in the top layer) to keep correct rho stratification for vertical eddy ! viscosity/diffusivity parameterization in PP82GLO. Later, to reduce ! pressure gradient errors associated with big mean density in underlying ! layers, we subtract the layer min rho from each layer. do k=1,k1 do j=2,j1 do i=2,i1 scr(i,j,k)=rho(i,j,k)-rhomin(k) enddo enddo enddo ! Use deviation of density from layer min density to get hydrostatic ! pressure gradient (reduces roundoff error and does not affect ! horizontal pressure gradient in the primitive hydroatic equations) ! Non-centered top half layer approximation: a good approximation ! because D(RHO)/DZ is very small in mixed layer ! tmpd=g*z(2) ! do j=2,j1 ! do i=2,i1 ! p(i,j,1)=p0(i,j)+tmpd*rho(i,j,1) ! enddo ! enddo ! do k=1,k2 ! tmpd=.5*g/odzw(k+1) ! do j=2,j1 ! do i=2,i1 ! p(i,j,k+1)=p(i,j,k)+tmpd*(rho(i,j,k)+rho(i,j,k+1)) ! enddo ! enddo ! enddo !!! 4th-order accurate 4th-order accurate do j=2,j1; jj=j-1; !(j:CELL, jj:SIZE OF WFACE, EASY TO READ) do i=2,i1; ii=i-1; !(i:CELL, ii:SIZE OF WFACE, EASY TO READ) wface(ii,jj,1)=p0(i,j) do k=1,k1 !(k: LAYER, TOP FACE: k, BOTTOM FACE: k+1) wface(ii,jj,k+1)=wface(ii,jj,k)+rho(i,j,k)*g/odz(k) end do end do end do ! C. CALCULATE VOLUME AVERAGED P (4TH-ORDER ACCURATE) do j=2,j1; jj=j-1; !(j:CELL, jj:SIZE OF WFACE, EASY TO READ) do i=2,i1; ii=i-1; !(i:CELL, ii:SIZE OF WFACE, EASY TO READ) do k=2,k2 !(k:INTERIOR LAYERS, TOP FACE: k, BOTTOM FACE: k+1) p(i,j,k)=(wface(ii,jj,k)+wface(ii,jj,k+1))*12.+ & (wface(ii,jj,k)+wface(ii,jj,k+1)-wface(ii,jj,k-1)-wface(ii,jj,k+2)) p(i,j,k)=p(i,j,k)*o24 end do ! 2ND-ORDER-ACCURATE IN THE TOP & BOTTOM LAYERS p(i,j,1) =.5*(wface(ii,jj,1) +wface(ii,jj,2)) p(i,j,k1)=.5*(wface(ii,jj,k1)+wface(ii,jj,k0)) end do do k=1,k1 p(1,j,k) =p(i1,j,k) p(i0,j,k)=p(2,j,k) end do end do end subroutine subroutine compute_pressure_gradient(k) ! first pressure differences !do j=2,j1 ! scratch array for 4th-order pressure gradient calculation ! avoid misusing data from overlying levels ! (put 0's where they are needed in xp) !do i=1,i1 ! xp(i,j)=iu(i,j,k)*(p(i+1,j,k)-p(i,j,k)) !enddo !enddo xp(1:i1,2:j1) = iu(1:i1,2:j1,k) * (p(2:i1+1,2:j1,k)-p(1:i1,2:j1,k)) call mpi_sendrecv(xp(2,1),1,m_evlon2d,m_w,1,xp(i0,1),1,m_evlon2d,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(xp(i2,1),1,m_evlon2d,m_e,1,xp(0,1),1,m_evlon2d,m_w,1,m_cart,istat,ierr) !do j=js,jf !do i=2,i1 ! yp(i,j)=iv(i,j,k)*(p(i,j+1,k)-p(i,j,k)) !enddo !enddo yp(2:i1,js:jf) = iv(2:i1,js:jf,k)*(p(2:i1,js+1:jf+1,k) - p(2:i1,js:jf,k)) call mpi_sendrecv(yp(1,2),1,m_evlat2d,m_s,1,yp(1,j0),1,m_evlat2d,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(yp(1,j2),1,m_evlat2d,m_n,1,yp(1,0),1,m_evlat2d,m_s,1,m_cart,istat,ierr) if (k .eq. 1) then do j=2,j1 do i=0,i0 xp(i,j)=iu0(i,j)*xp(i,j) !TS_NEW u(i,j,1)=iu0(i,j)*u(i,j,1) enddo enddo do j=js-1,jf+1 do i=2,i1 yp(i,j)=iv0(i,j)*yp(i,j) !TS_NEW v(i,j,1)=iv0(i,j)*v(i,j,1) enddo enddo endif ! basi! second order !do j=2,j1 !do i=2,i1 ! px(i,j)=6.*(xp(i,j)+xp(i-1,j)) ! py(i,j)=6.*(yp(i,j)+yp(i,j-1)) !enddo !enddo px(2:i1,2:j1) = 6.*(xp(2:i1,2:j1) + xp(1:i1-1,2:j1)) py(2:i1,2:j1) = 6.*(yp(2:i1,2:j1) + yp(2:i1,1:j1-1)) ! fourth order correction !do j=2,j1 ! do i=2,i1 ! px(i,j)=px(i,j)+iu(i-1,j,k)*iu(i,j,k)*(xp(i,j)+xp(i-1,j)-xp(i+1,j)-xp(i-2,j)) ! enddo !enddo px(2:i1,2:j1) = px(2:i1,2:j1)+iu(1:i1-1,2:j1,k)*iu(2:i1,2:j1,k)*(xp(2:i1,2:j1)+xp(1:i1-1,2:j1)-xp(3:i1+1,2:j1)-xp(0:i1-2,2:j1)) !do j=js+1,jf !do i=2,i1 ! py(i,j)=py(i,j)+iv(i,j-1,k)*iv(i,j,k)*(yp(i,j)+yp(i,j-1)-yp(i,j+1)-yp(i,j-2)) !enddo !enddo py(2:i1,js+1:jf)=py(2:i1,js+1:jf)+iv(2:i1,js:jf-1,k)*iv(2:i1,js+1:jf,k) & *(yp(2:i1,js+1:jf)+yp(2:i1,js:jf-1)-yp(2:i1,js+2:jf+1)-yp(2:i1,js-1:jf-2)) do j=2,j1 do i=2,i1 px(i,j)=o12*px(i,j) py(i,j)=o12*py(i,j) enddo enddo !px(2:i1, 2:j1) = o12 * px(2:i1, 2:j1) !py(2:i1, 2:j1) = o12 * py(2:i1, 2:j1) end subroutine subroutine compute_vertical_fluxes(k,lt,vdc) real*8, dimension(i0,j0,0:k0,2), intent(in) :: vdc ! --------------- ! vertical fluxes ! --------------- ! these are 4th-order-accurate in logical space ("stretched z-coord") ! which is logical when the stretching is near optimum. l=k+1 if (k==k1) then ! set bottom fluxes to zero (drag law applied after label 500) do j=1,j2 do i=1,i2 uz(i,j,lt)=0. vz(i,j,lt)=0. sz(i,j,lt)=0. tz(i,j,lt)=0. cfcz(i,j,lt)=0. enddo enddo !uz(1:i2, 1:j2, lt) = 0. !vz(1:i2, 1:j2, lt) = 0. !sz(1:i2, 1:j2, lt) = 0. !tz(1:i2, 1:j2, lt) = 0. !cfcz(1:i2, 1:j2, lt) = 0. else !do j=2,j1 ! do i=2,i1 ! scr(i,j,1)=6.*(u2(i,j,k)+u2(i,j,l)) ! scr(i,j,2)=6.*(v2(i,j,k)+v2(i,j,l)) ! scr(i,j,3)=6.*(s2(i,j,k)+s2(i,j,l)) ! scr(i,j,4)=6.*(t2(i,j,k)+t2(i,j,l)) ! scr(i,j,13)=6.*(c2(i,j,k)+c2(i,j,l)) ! enddo !enddo scr(2:i1, 2:j1, 1) = 6. * (u2(2:i1, 2:j1,k) + u2(2:i1, 2:j1,l)) scr(2:i1, 2:j1, 2) = 6. * (v2(2:i1, 2:j1,k) + v2(2:i1, 2:j1,l)) scr(2:i1, 2:j1, 3) = 6. * (s2(2:i1, 2:j1,k) + s2(2:i1, 2:j1,l)) scr(2:i1, 2:j1, 4) = 6. * (t2(2:i1, 2:j1,k) + t2(2:i1, 2:j1,l)) scr(2:i1, 2:j1,13) = 6. * (c2(2:i1, 2:j1,k) + c2(2:i1, 2:j1,l)) if (k/=1.and.k/=k2) then do j=2,j1 do i=2,i1 tmp=iw(i,j,k)*iw(i,j,k+2) scr(i,j,1)=scr(i,j,1)+tmp*(-ulf(i,j,k-1)+u2(i,j,k)+u2(i,j,l)-u2(i,j,k+2)) scr(i,j,2)=scr(i,j,2)+tmp*(-vlf(i,j,k-1)+v2(i,j,k)+v2(i,j,l)-v2(i,j,k+2)) scr(i,j,3)=scr(i,j,3)+tmp*(-slf(i,j,k-1)+s2(i,j,k)+s2(i,j,l)-s2(i,j,k+2)) scr(i,j,4)=scr(i,j,4)+tmp*(-tlf(i,j,k-1)+t2(i,j,k)+t2(i,j,l)-t2(i,j,k+2)) scr(i,j,13)=scr(i,j,13)+tmp*(-clf(i,j,k-1)+c2(i,j,k)+c2(i,j,l)-c2(i,j,k+2)) enddo enddo endif do n=1,4 do j=2,j1 !$omp do simd do i=2,i1 scr(i,j,n)=o12*scr(i,j,n) enddo enddo enddo do j=2,j1 !$omp do simd do i=2,i1 scr(i,j,13)=o12*scr(i,j,13) enddo enddo ! scr(2:i1,2:j1,1:4) = o12*scr(2:i1,2:j1,1:4) ! -1,j-1,lt)scr(2:i1,2:j1,13) = o12*scr(2:i1,2:j1,13) do j=2,j1 do i=2,i1 uz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,1)-ev(i-1,j-1,k)*(u1(i,j,l)-u1(i,j,k))*iw(i,j,l) vz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,2)-ev(i-1,j-1,k)*(v1(i,j,l)-v1(i,j,k))*iw(i,j,l) ! sz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,3)-0.1*hv(i-1,j-1,k)*(s1(i,j,l)-s1(i,j,k))*iw(i,j,l) ! tz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,4)-hv(i-1,j-1,k)*(t1(i,j,l)-t1(i,j,k))*iw(i,j,l) sz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,3)-vdc(i,j,k,2)*odzw(l)*(s1(i,j,l)-s1(i,j,k))*iw(i,j,l) tz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,4)-vdc(i,j,k,1)*odzw(l)*(t1(i,j,l)-t1(i,j,k))*iw(i,j,l) cfcz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,13)-hv(i-1,j-1,k)*(c1(i,j,l)-c1(i,j,k))*iw(i,j,l) enddo enddo endif ! test=0.d0 ! test=w(5,21,l)*scr(5,21,4)*iw(5,21,l)-vdc(5,21,k,1)*odzw(l)*(t1(5,21,l)-t1(5,21,k))*iw(5,21,l) ! if(myid.eq.0) write(99,*)'iw_k_l_tz_lt',iw(5,21,k),k,l,tz(4,20,lt),tz(4,20,lb) ! here ! if(myid.eq.0) write(99,*)'scr_w_vdc',scr(5,21,l),w(5,21,l),vdc(5,21,k,1) ! if(myid.eq.0) write(99,*)'test',test if (ldiag_budgets2) then !---------------------------- !20200615 output diagnostic !---------------------------- !WTS:Salt Flux Across Top Face, [gram/kilogram/s] !WTT:Heat Flux Across Top Face, [degC/s] if (k==1) then WTS(2:i1, 2:j1, 1) = sz(1:i2, 1:j2, 1)*odz(k) WTT(2:i1, 2:j1, 1) = tz(1:i2, 1:j2, 1)*odz(k) endif if (k==k1) then WTS(2:i1, 2:j1, l) = 0.0d0 WTT(2:i1, 2:j1, l) = 0.0d0 else do j=2,j1 do i=2,i1 WTS(i,j,l) = w(i,j,l)*scr(i,j,3)*odz(k) WTT(i,j,l) = w(i,j,l)*scr(i,j,4)*odz(k) enddo enddo endif endif !uz(1:i1-1,1:j1-1,lt)=w(2:i1,2:j1,l)*scr(2:i1,2:j1,1)-ev(1:i1-1,1:j1-1,k)*(u1(2:i1,2:j1,l)-u1(2:i1,2:j1,k))*iw(2:i1,2:j1,l) !vz(1:i1-1,1:j1-1,lt)=w(2:i1,2:j1,l)*scr(2:i1,2:j1,2)-ev(1:i1-1,1:j1-1,k)*(v1(2:i1,2:j1,l)-v1(2:i1,2:j1,k))*iw(2:i1,2:j1,l) !tz(1:i1-1,1:j1-1,lt)=w(2:i1,2:j1,l)*scr(2:i1,2:j1,4)-hv(1:i1-1,1:j1-1,k)*(t1(2:i1,2:j1,l)-t1(2:i1,2:j1,k))*iw(2:i1,2:j1,l) !cfcz(1:i1-1,1:j1-1,lt)=w(2:i1,2:j1,l)*scr(2:i1,2:j1,13)-hv(1:i1-1,1:j1-1,k)*(c1(2:i1,2:j1,l)-c1(2:i1,2:j1,k))*iw(2:i1,2:j1,l) ! set scr(3,j,5) to store u2(3,j,k) so as v2,s2 and t2 call mpi_sendrecv(u2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,5),1,m_vlon2d,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(v2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,6),1,m_vlon2d,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(s2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,7),1,m_vlon2d,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(t2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,8),1,m_vlon2d,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(c2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,14),1,m_vlon2d,m_e,1,m_cart,istat,ierr) ! set scr(i2,j,5) to store u2(i2,j,k) so as v2,s2 and t2 call mpi_sendrecv(u2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,5),1,m_vlon2d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(v2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,6),1,m_vlon2d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(s2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,7),1,m_vlon2d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(t2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,8),1,m_vlon2d,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(c2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,14),1,m_vlon2d,m_w,1,m_cart,istat,ierr) ! set scr(i,3,9) to store u2(i,3,k) so as v2,s2 and t2 call mpi_sendrecv(u2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,9),1,m_vlat2d,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(v2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,10),1,m_vlat2d,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(s2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,11),1,m_vlat2d,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(t2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,12),1,m_vlat2d,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(c2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,15),1,m_vlat2d,m_n,1,m_cart,istat,ierr) ! set scr(i,j2,9) to store u2(i,j2,k) so as v2,s2 and t2 call mpi_sendrecv(u2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,9),1,m_vlat2d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(v2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,10),1,m_vlat2d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(s2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,11),1,m_vlat2d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(t2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,12),1,m_vlat2d,m_s,1,m_cart,istat,ierr) call mpi_sendrecv(c2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,15),1,m_vlat2d,m_s,1,m_cart,istat,ierr) end subroutine subroutine compute_longitudinal_fluxes(k) ! ------------------- ! longitudinal fluxes ! ------------------- do j=2,j1 ! higher order interpolations !do i=1,i1 ! scr(i,j,1)=6.*(u2(i,j,k)+u2(i+1,j,k)) ! scr(i,j,2)=6.*(v2(i,j,k)+v2(i+1,j,k)) ! scr(i,j,3)=6.*(s2(i,j,k)+s2(i+1,j,k)) ! scr(i,j,4)=6.*(t2(i,j,k)+t2(i+1,j,k)) ! scr(i,j,13)=6.*(c2(i,j,k)+c2(i+1,j,k)) !enddo scr(1:i1, j, 1) = 6. * (u2(1:i1, j, k) + u2(2:i1+1, j, k)) scr(1:i1, j, 2) = 6. * (v2(1:i1, j, k) + v2(2:i1+1, j, k)) scr(1:i1, j, 3) = 6. * (s2(1:i1, j, k) + s2(2:i1+1, j, k)) scr(1:i1, j, 4) = 6. * (t2(1:i1, j, k) + t2(2:i1+1, j, k)) scr(1:i1, j,13) = 6. * (c2(1:i1, j, k) + c2(2:i1+1, j, k)) tmp=iu(0,j,k)*iu(2,j,k) scr(1,j,1)=scr(1,j,1)+tmp*(-scr(i2,j,5)+u2(1,j,k)+u2(2,j,k)-u2(3,j,k)) scr(1,j,2)=scr(1,j,2)+tmp*(-scr(i2,j,6)+v2(1,j,k)+v2(2,j,k)-v2(3,j,k)) scr(1,j,3)=scr(1,j,3)+tmp*(-scr(i2,j,7)+s2(1,j,k)+s2(2,j,k)-s2(3,j,k)) scr(1,j,4)=scr(1,j,4)+tmp* (-scr(i2,j,8)+t2(1,j,k)+t2(2,j,k)-t2(3,j,k)) scr(1,j,13)=scr(1,j,13)+tmp*(-scr(i2,j,14)+c2(1,j,k)+c2(2,j,k)-c2(3,j,k)) tmp=iu(i2,j,k)*iu(i0,j,k) scr(i1,j,1)=scr(i1,j,1)+tmp*(-u2(i2,j,k)+u2(i1,j,k)+u2(i0,j,k)-scr(3,j,5)) scr(i1,j,2)=scr(i1,j,2)+tmp*(-v2(i2,j,k)+v2(i1,j,k)+v2(i0,j,k)-scr(3,j,6)) scr(i1,j,3)=scr(i1,j,3)+tmp*(-s2(i2,j,k)+s2(i1,j,k)+s2(i0,j,k)-scr(3,j,7)) scr(i1,j,4)=scr(i1,j,4)+tmp*(-t2(i2,j,k)+t2(i1,j,k)+t2(i0,j,k)-scr(3,j,8)) scr(i1,j,13)=scr(i1,j,13)+tmp*(-c2(i2,j,k)+c2(i1,j,k)+c2(i0,j,k)-scr(3,j,14)) do i=2,i2 tmp=iu(i-1,j,k)*iu(i+1,j,k) scr(i,j,1)=scr(i,j,1)+tmp*(-u2(i-1,j,k)+u2(i,j,k)+u2(i+1,j,k)-u2(i+2,j,k)) scr(i,j,2)=scr(i,j,2)+tmp*(-v2(i-1,j,k)+v2(i,j,k)+v2(i+1,j,k)-v2(i+2,j,k)) scr(i,j,3)=scr(i,j,3)+tmp*(-s2(i-1,j,k)+s2(i,j,k)+s2(i+1,j,k)-s2(i+2,j,k)) scr(i,j,4)=scr(i,j,4)+tmp*(-t2(i-1,j,k)+t2(i,j,k)+t2(i+1,j,k)-t2(i+2,j,k)) scr(i,j,13)=scr(i,j,13)+tmp*(-c2(i-1,j,k)+c2(i,j,k)+c2(i+1,j,k)-c2(i+2,j,k)) enddo do n=1,4 !$omp do simd do i=1,i1 scr(i,j,n)=o12*scr(i,j,n) enddo enddo !$omp do simd do i=1,i1 scr(i,j,13)=o12*scr(i,j,13) enddo ! scr(1:i1, j, 1:4) = o12 * scr(1:i1, j, 1:4) ! scr(1:i1, j, 13) = o12 * scr(1:i1, j, 13) do i=1,i1 am=dmx(i,j,k)*odx(j) ah=dhx(i,j)*odx(j) ! ux could be simply u(i,j,k)**2 + ...., but this does not conserve ! energy and is unstable ux(i,j-1)=(u(i,j,k)*scr(i,j,1)-am*(u1(i+1,j,k)-u1(i,j,k)))*iu(i,j,k) vx(i,j-1)=(u(i,j,k)*scr(i,j,2)-am*(v1(i+1,j,k)-v1(i,j,k)))*iu(i,j,k) sx(i,j-1)=(u(i,j,k)*scr(i,j,3)-ah*(s1(i+1,j,k)-s1(i,j,k)))*iu(i,j,k) tx(i,j-1)=(u(i,j,k)*scr(i,j,4)-ah*(t1(i+1,j,k)-t1(i,j,k)))*iu(i,j,k) cfcx(i,j-1)=(u(i,j,k)*scr(i,j,13)-ah*(c1(i+1,j,k)-c1(i,j,k)))*iu(i,j,k) if (ldiag_budgets2) then !-------------------------------------- !20200615 output diagnostic:HDIFT,HDIFS !-------------------------------------- sx_diffu(i,j-1)=(-ah*(s1(i+1,j,k)-s1(i,j,k)))*iu(i,j,k) tx_diffu(i,j-1)=(-ah*(t1(i+1,j,k)-t1(i,j,k)))*iu(i,j,k) endif enddo enddo if (ldiag_budgets2) then !---------------------------- !20200615 output diagnostic !---------------------------- !UES:Salt Flux in grid-x direction, [gram/kilogram/s] !UET:Heat Flux in grid-x direction, [degC/s] do j=2,j1 do i=1,i1 UES(i,j,k) = u(i,j,k)*scr(i,j,3)*iu(i,j,k)*odx(j) UET(i,j,k) = u(i,j,k)*scr(i,j,4)*iu(i,j,k)*odx(j) enddo enddo endif end subroutine subroutine compute_latitudinal_fluxes(k) !$OMP DECLARE SIMD (compute_latitudinal_fluxes) UNIFORM(k) ! ------------------ ! latitudinal fluxes ! ------------------ !do j=1,jf ! do i=2,i1 ! scr(i,j,1)=6.*(u2(i,j,k)+u2(i,j+1,k)) ! scr(i,j,2)=6.*(v2(i,j,k)+v2(i,j+1,k)) ! scr(i,j,3)=6.*(s2(i,j,k)+s2(i,j+1,k)) ! scr(i,j,4)=6.*(t2(i,j,k)+t2(i,j+1,k)) ! scr(i,j,13)=6.*(c2(i,j,k)+c2(i,j+1,k)) ! enddo !enddo scr(2:i1, 1:jf, 1) = 6. * (u2(2:i1, 1:jf, k) + u2(2:i1, 2:jf+1, k)) scr(2:i1, 1:jf, 2) = 6. * (v2(2:i1, 1:jf, k) + v2(2:i1, 2:jf+1, k)) scr(2:i1, 1:jf, 3) = 6. * (s2(2:i1, 1:jf, k) + s2(2:i1, 2:jf+1, k)) scr(2:i1, 1:jf, 4) = 6. * (t2(2:i1, 1:jf, k) + t2(2:i1, 2:jf+1, k)) scr(2:i1, 1:jf,13) = 6. * (c2(2:i1, 1:jf, k) + c2(2:i1, 2:jf+1, k)) do j=2,j2 do i=2,i1 tmp=iv(i,j-1,k)*iv(i,j+1,k) scr(i,j,1)=scr(i,j,1)+tmp*(-u2(i,j-1,k)+u2(i,j,k)+u2(i,j+1,k)-u2(i,j+2,k)) scr(i,j,2)=scr(i,j,2)+tmp*(-v2(i,j-1,k)+v2(i,j,k)+v2(i,j+1,k)-v2(i,j+2,k)) scr(i,j,3)=scr(i,j,3)+tmp*(-s2(i,j-1,k)+s2(i,j,k)+s2(i,j+1,k)-s2(i,j+2,k)) scr(i,j,4)=scr(i,j,4)+tmp*(-t2(i,j-1,k)+t2(i,j,k)+t2(i,j+1,k)-t2(i,j+2,k)) scr(i,j,13)=scr(i,j,13)+tmp*(-c2(i,j-1,k)+c2(i,j,k)+c2(i,j+1,k)-c2(i,j+2,k)) enddo enddo do i=2,i1 tmp=iv(i,0,k)*iv(i,2,k) scr(i,1,1)=scr(i,1,1)+tmp*(-scr(i,j2,9)+u2(i,1,k)+u2(i,2,k)-u2(i,3,k)) scr(i,1,2)=scr(i,1,2)+tmp*(-scr(i,j2,10)+v2(i,1,k)+v2(i,2,k)-v2(i,3,k)) scr(i,1,3)=scr(i,1,3)+tmp*(-scr(i,j2,11)+s2(i,1,k)+s2(i,2,k)-s2(i,3,k)) scr(i,1,4)=scr(i,1,4)+tmp*(-scr(i,j2,12)+t2(i,1,k)+t2(i,2,k)-t2(i,3,k)) scr(i,1,13)=scr(i,1,13)+tmp*(-scr(i,j2,15)+c2(i,1,k)+c2(i,2,k)-c2(i,3,k)) enddo do i=2,i1 tmp=iv(i,j2,k)*iv(i,j0,k) scr(i,j1,1)=scr(i,j1,1)+tmp*(-u2(i,j2,k)+u2(i,j1,k)+u2(i,j0,k)-scr(i,3,9)) scr(i,j1,2)=scr(i,j1,2)+tmp*(-v2(i,j2,k)+v2(i,j1,k)+v2(i,j0,k)-scr(i,3,10)) scr(i,j1,3)=scr(i,j1,3)+tmp*(-s2(i,j2,k)+s2(i,j1,k)+s2(i,j0,k)-scr(i,3,11)) scr(i,j1,4)=scr(i,j1,4)+tmp*(-t2(i,j2,k)+t2(i,j1,k)+t2(i,j0,k)-scr(i,3,12)) scr(i,j1,13)=scr(i,j1,13)+tmp*(-c2(i,j2,k)+c2(i,j1,k)+c2(i,j0,k)-scr(i,3,15)) enddo do n=1,4 do j=js,jf !$omp do simd do i=2,i1 scr(i,j,n)=o12*scr(i,j,n) enddo enddo enddo do j=js,jf !$omp do simd do i=2,i1 scr(i,j,13)=o12*scr(i,j,13) enddo enddo !do j=js,jf ! scr(2:i1, j, 13) = o12 * scr(2:i1, j, 13) !enddo !scr(2:i1, js:jf, 1:4) = o12 * scr(2:i1, js:jf, 1:4) !scr(2:i1, js:jf, 13) = o12 * scr(2:i1, js:jf, 13) do j=js,jf do i=2,i1 am=dmy(i,j,k)*odyv(j) ah=dhy(i,j)*odyv(j) uy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,1)-am*(u1(i,j+1,k)-u1(i,j,k)))*iv(i,j,k) vy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,2)-am*(v1(i,j+1,k)-v1(i,j,k)))*iv(i,j,k) sy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,3)-ah*(s1(i,j+1,k)-s1(i,j,k)))*iv(i,j,k) ty(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,4)-ah*(t1(i,j+1,k)-t1(i,j,k)))*iv(i,j,k) cfcy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,13)-ah*(c1(i,j+1,k)-c1(i,j,k)))*iv(i,j,k) if (ldiag_budgets2) then !-------------------------------------- !20200615 output diagnostic:HDIFT,HDIFS !-------------------------------------- sy_diffu(i-1,j)=csv(j)*(-ah*(s1(i,j+1,k)-s1(i,j,k)))*iv(i,j,k) ty_diffu(i-1,j)=csv(j)*(-ah*(t1(i,j+1,k)-t1(i,j,k)))*iv(i,j,k) endif enddo enddo if (ldiag_budgets2) then !---------------------------- !20200615 output diagnostic !---------------------------- !VNS:Salt Flux in grid-y direction, [gram/kilogram/s] !VNT:Heat Flux in grid-y direction, [degC/s] do j=js,jf do i=2,i1 odyj=ocs(j)*ody(j) VNS(i,j,k) = csv(j)*v(i,j,k)*scr(i,j,3)*iv(i,j,k)*odyj VNT(i,j,k) = csv(j)*v(i,j,k)*scr(i,j,4)*iv(i,j,k)*odyj enddo enddo endif end subroutine subroutine init_prepexec !exchanging the boundary points from inner domain prep.data (i2,j2) !! WORKING !! !rundata !zkb !init_solver (at "i2+1" only) !set_seasonal_data are exchenged within its own subroutines end subroutine end module