subroutine barotp_dry_and_nodry(m,n) use mod_xc ! HYCOM communication interface use mod_pipe ! HYCOM debugging interface use mod_cb_arrays ! HYCOM saved arrays c c --- ------------------------------------------------------------------------ c --- advance barotropic equations from baroclinic time level -m- to level -n c --- ------------------------------------------------------------------------ c --- ------------------------------------------------------------------------ c --- advance barotropic equations. c --- on entry: -n- is time t-dt, -m- is time t c --- on exit: -m- is time t, -n- is time t+dt c --- time level 3 is only used internally (n and m are always 1 or 2). c c --- LeapFrog version based on: c --- Y. Morel, Baraille, R., Pichon A. (2008) "Time splitting and c --- linear stability of the slow part of the barotropic component", c --- Ocean Modeling, 23, pp 73-81. c c --- Wetting and Drying only at needed points (where depths < depth_wet) c c===========================================================--============ c implicit none c real z1, wblpf c integer m, n, ml, nl, mn integer lstep1, icof integer i, j, l, k, lll, ia, ib, ja, jb,margin c c Mass conservation #if defined(RELO) real, save, allocatable, dimension(:,:) :: & pbavo,ubavo,vbavo, & flxloc,flyloc,uflxba,vflxba c if (.not.allocated(pbavo)) then allocate( & pbavo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & ubavo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vbavo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & flxloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & flyloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & uflxba(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vflxba(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) call mem_stat_add( 6*(idm+2*nbdy)*(jdm+2*nbdy) ) pbavo = r_init ubavo = r_init vbavo = r_init flxloc = r_init flyloc = r_init uflxba = r_init vflxba = r_init endif #else real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & pbavo,ubavo,vbavo, & flxloc,flyloc,uflxba,vflxba #endif c c======================================================================= c if (btrlfr .and. nstep .gt. 1) then lstep1 = lstep+lstep icof = 2 else lstep1 = lstep icof = 1 endif c c --- utotn,vtotn from momtum is time step t-1 to t+1 barotropic tendency call xctilr(utotn( 1-nbdy,1-nbdy ),1, 1, 6,6, halo_uv) call xctilr(vtotn( 1-nbdy,1-nbdy ),1, 1, 6,6, halo_vv) c margin = nbdy ! do j=1,jj ! do i=1,ii ! if (SEA_P) then ! onetao(i,j,n) = oneta(i,j,n) ! endif ! enddo ! enddo ! uflxba(:,:) = 0. ! vflxba(:,:) = 0. cCL/RB MPI bug correction 2011-01 !!Alex call xctilr(oneta( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_ps) call xctilr(oneta( 1-nbdy,1-nbdy, m),1, 1, 6,6, halo_ps) if (btrmas) then uflxba(:,:) = 0.0 vflxba(:,:) = 0.0 flxloc(:,:) = 0.0 flyloc(:,:) = 0.0 onetao(:,:,n) = oneta(:,:,n) endif c c Save t-1 for RA filter c ---------------------- !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) ubavo(i,j) = (1.-ra2fac)*ubavg(i,j,m)+.5*ra2fac*ubavg(i,j,n) enddo enddo do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vbavo(i,j) = (1.-ra2fac)*vbavg(i,j,m)+.5*ra2fac*vbavg(i,j,n) enddo enddo do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pbavo(i,j) = (1.-ra2fac)*pbavg(i,j,m)+.5*ra2fac*pbavg(i,j,n) enddo enddo enddo c if (btrlfr) then !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) ubavg(i,j,m) = ubavg(i,j,n) enddo enddo do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vbavg(i,j,m) = vbavg(i,j,n) enddo enddo do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pbavg(i,j,m) = pbavg(i,j,n) enddo enddo enddo c wbaro is set to zero for the first time step c -------------------------------------------- wblpf = 0. else wblpf = wbaro endif c ml = n : t-1, mn = 3 : t, nl = m : t+1 c - - - - - - - - - - - - - - - - - - - ml = n mn = 3 nl = m c do 840 lll = 1, lstep1 c i = ml ml = mn mn = nl nl = i c call xctilr(pbavg( 1-nbdy,1-nbdy,mn ),1, 1, 4,4, halo_ps) call xctilr(ubavg( 1-nbdy,1-nbdy,mn ),1, 1, 3,3, halo_uv) call xctilr(vbavg( 1-nbdy,1-nbdy,mn ),1, 1, 3,3, halo_vv) call xctilr(ubavg( 1-nbdy,1-nbdy,ml ),1, 1, 1,1, halo_uv) call xctilr(vbavg( 1-nbdy,1-nbdy,ml ),1, 1, 1,1, halo_vv) c c Continuity c ========== c Boundary conditions c -------------------- #if ! defined(RELO) if (lbflag.eq.1) then call latbdp( nl) elseif (lbflag.eq.3) then call latbdf( nl,lll+1) endif #endif if (lbflag.eq.2) then call latbdt( nl,lll) elseif (lbflag.eq.4) then call latbdtf(nl,lll) endif c margin = 3 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) util3(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,mn)) enddo enddo enddo c first order u-flux for the wet./dry. area c ----------------------------------------- margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) if(ubavg(i,j,mn).ge.0.) then uflux(i,j) = util3(i-1,j)*ubavg(i,j,mn)*scuy(i,j) else uflux(i,j) = util3(i,j)*ubavg(i,j,mn)*scuy(i,j) endif uflux2(i,j) = 0.5*ubavg(i,j,mn)*scuy(i,j)* . (util3(i,j)+util3(i-1,j)) - uflux(i,j) flxloc(i,j) = uflux(i,j) enddo enddo enddo c Second order u-flux for the non wet./dry. area c ---------------------------------------------- margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_nw(j) do i=max(1-margin,ifu_nw(j,l)),min(ii+margin,ilu_nw(j,l)) uflux(i,j) = 0.5*ubavg(i,j,mn)*scuy(i,j)* . (util3(i,j)+util3(i-1,j)) flxloc(i,j) = uflux(i,j) enddo enddo enddo c Boundary conditions c ------------------- margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) i=max(1-nbdy,ifu(j,l)-1) if (iuopn(i,j).ne.0) then c In the SHOM version, boundary values are stored in pbsavu, pbsavv c Probably stored in pbavg(.,mn) in the standard version c For SHOM uflux(i,j)=max(0.,pbot(i,j)*onetai(i,j)+pbsavu(i,j)) uflux(i,j)=max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,mn)) . *ubavg(i,j,mn)*scuy(i,j) else uflux(i,j)=0. endif flxloc(i,j) = uflux(i,j) uflux2(i,j) = 0. i=min(ii+nbdy,ilu(j,l)+1) if (iuopn(i,j).ne.0) then c For SHOM uflux(i,j)=max(0.,pbot(i-1,j)*onetai(i-1,j)+pbsavu(i,j)) uflux(i,j)=max(0.,pbot(i-1,j)*onetai(i-1,j)+ . pbavg(i,j,mn)) . *ubavg(i,j,mn)*scuy(i,j) else uflux(i,j)=0. endif flxloc(i,j) = uflux(i,j) uflux2(i,j) = 0. enddo enddo c first order v-flux for the wet./dry. area c ----------------------------------------- !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) if(vbavg(i,j,mn).ge.0.) then vflux(i,j)=util3(i,j-1)*vbavg(i,j,mn)*scvx(i,j) else vflux(i,j)=util3(i,j)*vbavg(i,j,mn)*scvx(i,j) endif vflux2(i,j) = 0.5*scvx(i,j)*vbavg(i,j,mn)* . (util3(i,j)+util3(i,j-1)) - vflux(i,j) flyloc(i,j) = vflux(i,j) enddo enddo enddo c Second order v-flux for the non wet./dry. area c ---------------------------------------------- margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_nw(j) do i=max(1-margin,ifv_nw(j,l)),min(ii+margin,ilv_nw(j,l)) vflux(i,j) = 0.5*scvx(i,j)*vbavg(i,j,mn)* . (util3(i,j)+util3(i,j-1)) flyloc(i,j) = vflux(i,j) enddo enddo enddo c Boundary conditions c ------------------- margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) do i=1-margin,ii+margin do l=1,jsv(i) j=max(1-nbdy,jfv(i,l)-1) if (ivopn(i,j).ne.0) then ccc For SHOM vflux(i,j)=max(0.,pbot(i,j)*onetai(i,j)+pbsavv(i,j)) vflux(i,j)=max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,mn)) . *vbavg(i,j,mn)*scvx(i,j) else vflux(i,j)=0. endif flyloc(i,j) = vflux(i,j) vflux2(i,j)=0. j=min(jj+nbdy,jlv(i,l)+1) if (ivopn(i,j).ne.0) then ccc For SHOM vflux(i,j)=max(0.,pbot(i,j-1)*onetai(i,j-1)+pbsavv(i,j)) vflux(i,j)=max(0.,pbot(i,j-1)*onetai(i,j-1) + . pbavg(i,j,mn)) . *vbavg(i,j,mn)*scvx(i,j) else vflux(i,j)=0. endif flyloc(i,j) = vflux(i,j) vflux2(i,j)=0. enddo enddo c Integration c ----------- margin = 0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pbavg(i,j,nl)= . ((1.-wblpf)*pbavg(i,j,mn)+ . wblpf *pbavg(i,j,ml) ) . - (uflux(i+1,j)-uflux(i,j) . + vflux(i,j+1)-vflux(i,j) . )*dlt*scp2i(i,j)*(1.+wblpf) enddo enddo enddo call xctilr(pbavg( 1-nbdy,1-nbdy,nl ),1, 1, 3,3, halo_ps) c c End of the FCT scheme for the wet./dry. area c -------------------------------------------- margin = 3 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp_w(j) do i=max(1-margin,ifp_w(j,l)),min(ii+margin,ilp_w(j,l)) util4(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) enddo i = max(1-nbdy,ifp_w(j,l)-1) if (ip(i,j).eq.1) . util4(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) i = min(ii+nbdy,ilp_w(j,l)+1) if (ip(i,j).eq.1) . util4(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) enddo enddo do i = 1-margin, ii+margin do l = 1, jsp_w(i) j = max(1-nbdy,jfp_w(i,l)-1) if (ip(i,j).eq.1) . util4(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) j = min(jj+nbdy,jlp_w(i,l)+1) if (ip(i,j).eq.1) . util4(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) enddo enddo c R+,R- calculations c ------------------ margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i,ia,ib,ja,jb) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp_w(j) do i=max(1-margin,ifp_w(j,l)),min(ii+margin,ilp_w(j,l)) ia=i-1 if (ip(ia,j).eq.0) ia=i ib=i+1 if (ip(ib,j).eq.0) ib=i ja=j-1 if (ip(i,ja).eq.0) ja=j jb=j+1 if (ip(i,jb).eq.0) jb=j util1(i,j)=max(util4(i,j),util4(ia,j),util4(ib,j), . util4(i,ja),util4(i,jb)) util2(i,j)=max(0., . min(util4(i,j),util4(ia,j),util4(ib,j), . util4(i,ja),util4(i,jb))) c util1(i,j)=(util1(i,j)-util4(i,j)) . /(((max(0.,uflux2(i,j))-min(0.,uflux2(i+1,j))) . +(max(0.,vflux2(i,j))-min(0.,vflux2(i,j+1)))+epsil) . *dlt*(1.+wblpf)*scp2i(i,j)) c util2(i,j)=(util2(i,j)-util4(i,j)) . /(((min(0.,uflux2(i,j))-max(0.,uflux2(i+1,j))) . +(min(0.,vflux2(i,j))-max(0.,vflux2(i,j+1)))-epsil) . *dlt*(1.+wblpf)*scp2i(i,j)) enddo enddo enddo c Final fluxes c ------------ margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin c do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) if (uflux2(i,j).ge.0.) then uflux(i,j) = uflux2(i,j) . *min(1.,util1(i,j),util2(i-1,j)) else uflux(i,j) = uflux2(i,j) . *min(1.,util2(i,j),util1(i-1,j)) endif flxloc(i,j) = flxloc(i,j) + uflux(i,j) enddo enddo c do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) if (vflux2(i,j).ge.0.) then vflux(i,j) = vflux2(i,j) . *min(1.,util1(i,j),util2(i,j-1)) else vflux(i,j) = vflux2(i,j) . *min(1.,util2(i,j),util1(i,j-1)) endif flyloc(i,j) = flyloc(i,j) + vflux(i,j) enddo enddo enddo c set fluxes to 0 at the interface c -------------------------------- !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) i=max(1-nbdy,ifu_w(j,l)-1) uflux(i,j)=0. i=min(ii+nbdy,ilu_w(j,l)+1) uflux(i,j)=0. enddo enddo c !$OMP PARALLEL DO PRIVATE(j,l,i) do i=1-margin,ii+margin do l=1,jsv_w(i) j=max(1-nbdy,jfv_w(i,l)-1) vflux(i,j)=0. j=min(jj+nbdy,jlv_w(i,l)+1) vflux(i,j)=0. enddo enddo c Second order integration for the wet./dry.area c ---------------------------------------------- margin = 0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp_w(j) do i=max(1-margin,ifp_w(j,l)),min(ii+margin,ilp_w(j,l)) pbavg(i,j,nl)=pbavg(i,j,nl) . - ((iu_w(i+1,j)*uflux(i+1,j)-iu_w(i,j)*uflux(i,j))+ . (iv_w(i,j+1)*vflux(i,j+1)-iv_w(i,j)*vflux(i,j))) . *dlt*(1.+wblpf)*scp2i(i,j) enddo enddo enddo c Barotropic fluxes for mass conservation c --------------------------------------- do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) uflxba(i,j) = uflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flxloc(i,j) enddo enddo enddo do j=1-margin,jj+margin do l=1,isu(j) i=max(1-nbdy,ifu(j,l)-1) uflxba(i,j) = uflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flxloc(i,j) i=min(ii+nbdy,ilu(j,l)+1) uflxba(i,j) = uflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flxloc(i,j) enddo enddo do j=1-margin,jj+margin do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vflxba(i,j) = vflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flyloc(i,j) enddo enddo enddo do i=1-margin,ii+margin do l=1,jsv(i) j=max(1-nbdy,jfv(i,l)-1) vflxba(i,j) = vflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flyloc(i,j) j=min(jj+nbdy,jlv(i,l)+1) vflxba(i,j) = vflxba(i,j) + coeflx(lll+1,icof)*(1.+wblpf) . * flyloc(i,j) enddo enddo c call xctilr(pbavg( 1-nbdy,1-nbdy,nl ),1, 1, 3,3, halo_ps) c c Momentum equations c ================== c u-equation c - - - - - - margin = 2 c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) util2(i,j) = (-utotn(i,j) . +thref*(pbavg(i,j,nl)-pbavg(i-1,j,nl))*scuxi(i,j) . -((vbavg(i,j,mn)*depthv(i,j) . +vbavg(i,j+1,mn)*depthv(i,j+1))+ . (vbavg(i-1,j,mn)*depthv(i-1,j) . +vbavg(i-1,j+1,mn)*depthv(i-1,j+1)))* . (0.125*(pvtrop(i,j)+pvtrop(i,j+1)))) enddo enddo enddo c Integration for non wet./wry.area c --------------------------------- margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_nw(j) do i=max(1-margin,ifu_nw(j,l)),min(ii+margin,ilu_nw(j,l)) ubavg(i,j,nl) = ((1.-wblpf)*ubavg(i,j,mn)+ . wblpf*ubavg(i,j,ml))-(1.+wblpf)*dlt*util2(i,j) enddo enddo enddo c margin = 3 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp_w(j) do i=max(1-margin,ifp_w(j,l)-1),min(ii+margin,ilp_w(j,l)+1) util3(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) enddo enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do i=1-margin,ii+margin do l=1,jsp_w(i) j=max(1-margin,jfp_w(i,l)-1) util3(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) j=min(jj+margin,jlp_w(i,l)+1) util3(i,j) = max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,nl)) enddo enddo c Weights c ------- margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)-1),min(ii+margin,ilu_w(j,l)+1) util1(i,j)=max(0., . min(util3(i,j),util3(i-1,j),hcpu(i,j))) util2(i,j)=util1(i,j)*util2(i,j) enddo enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do i=1-margin,ii+margin do l=1,jsu_w(i) j = max(1-margin,jfu_w(i,l)-1) if (iu_w(i,j).ne.1) then !!Alex add iu_w from Remy util1(i,j)=max(0.,min(util3(i,j),util3(i-1,j),hcpu(i,j))) util2(i,j)=util1(i,j)*util2(i,j) endif j = min(jj+margin,jlu_w(i,l)+1) if (iu_w(i,j).ne.1) then !!Alex add iu_w from Remy util1(i,j)=max(0.,min(util3(i,j),util3(i-1,j),hcpu(i,j))) util2(i,j)=util1(i,j)*util2(i,j) endif enddo enddo c c margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) ia=i-1 if (iu(ia,j).ne.1) ia=i ib=i+1 if (iu(ib,j).ne.1) ib=i ja=j-1 if (iu(i,ja).ne.1) ja=j jb=j+1 if (iu(i,jb).ne.1) jb=j util4(i,j)=(util2(i,j)+(hcpu(i,j)-util1(i,j))* . (util2 (ia,j)+util2 (ib,j)+util2 (i,ja)+util2 (i,jb))/ . (util1(ia,j)+util1(ib,j)+util1(i,ja)+util1(i,jb)+ . epspu(i,j)))*invhcpu(i,j) enddo enddo enddo c Integration for wet./wry.area c ----------------------------- util5(:,:) = 0. ! Probably not necessary - to be checked later c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) util5(i,j)= ((1.-wblpf)*ubavg(i,j,mn)+ . wblpf*ubavg(i,j,ml))-(1.+wblpf)*dlt*util4(i,j) enddo enddo enddo c v-equation c - - - - - - c margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) util2(i,j) = (-vtotn(i,j) . +thref*(pbavg(i,j,nl)-pbavg(i,j-1,nl))*scvyi(i,j) . +((ubavg(i,j,mn)*depthu(i,j) . +ubavg(i+1,j,mn)*depthu(i+1,j))+ . (ubavg(i,j-1,mn)*depthu(i,j-1) . +ubavg(i+1,j-1,mn)*depthu(i+1,j-1)))* . (0.125*(pvtrop(i,j)+pvtrop(i+1,j)))) enddo enddo enddo c Integration for non wet./wry.area c --------------------------------- margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_nw(j) do i=max(1-margin,ifv_nw(j,l)),min(ii+margin,ilv_nw(j,l)) vbavg(i,j,nl) = ((1.-wblpf)*vbavg(i,j,mn)+ . wblpf*vbavg(i,j,ml))-(1.+wblpf)*dlt*util2(i,j) enddo enddo enddo c Weights c ------- margin = 2 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)-1),min(ii+margin,ilv_w(j,l)+1) util1(i,j)=max(0., . min(util3(i,j-1),util3(i,j),hcpv(i,j))) util2(i,j)=util1(i,j)*util2(i,j) enddo enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do i=1-margin,ii+margin !!Alex bug corrected from jj to ii do l=1,jsv_w(i) j=max(1-margin,jfv_w(i,l)-1) if (iv_w(i,j).ne.1) then util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcpv(i,j))) util2(i,j)=util1(i,j)*util2(i,j) endif j=min(jj+margin,jlv_w(i,l)+1) if (iv_w(i,j).ne.1) then util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcpv(i,j))) util2(i,j)=util1(i,j)*util2(i,j) endif enddo enddo c margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) ia=i-1 if (iv(ia,j).ne.1) ia=i ib=i+1 if (iv(ib,j).ne.1) ib=i ja=j-1 if (iv(i,ja).ne.1) ja=j jb=j+1 if (iv(i,jb).ne.1) jb=j util4(i,j)=(util2(i,j)+(hcpv(i,j)-util1(i,j))* . (util2 (ia,j)+util2 (ib,j)+util2 (i,ja)+util2 (i,jb))/ . (util1(ia,j)+util1(ib,j)+util1(i,ja)+util1(i,jb)+ . epspv(i,j)))*invhcpv(i,j) enddo enddo enddo c Integration for wet./wry.area c ----------------------------- util6(:,:) = 0. ! Probably not necessary - to be checked later c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) util6(i,j)= ((1.-wblpf)*vbavg(i,j,mn)+ . wblpf*vbavg(i,j,ml))-(1.+wblpf)*dlt*util4(i,j) enddo enddo enddo c Velocity relaxation when dp is low c ================================== margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) util1(i,j)=min(hcu(i,j),util3(i,j),util3(i-1,j)) enddo i = max(1-margin,ifu_w(j,l)-1) util1(i,j)=min(hcu(i,j),util3(i,j),util3(i-1,j)) if (iu_nw(i,j).eq.1) util5(i,j)=ubavg(i,j,nl) i = min(ii+margin,ilu_w(j,l)+1) util1(i,j)=min(hcu(i,j),util3(i,j),util3(i-1,j)) if (iu_nw(i,j).eq.1) util5(i,j)=ubavg(i,j,nl) enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do i=1-margin,ii+margin do l=1,jsu_w(i) j = max(1-margin,jfu_w(i,l)-1) util1(i,j)=max(0.,min(util3(i,j),util3(i-1,j),hcu(i,j))) if (iu_nw(i,j).eq.1) util5(i,j)=ubavg(i,j,nl) j = min(jj+margin,jlu_w(i,l)+1) util1(i,j)=max(0.,min(util3(i,j),util3(i-1,j),hcu(i,j))) if (iu_nw(i,j).eq.1) util5(i,j)=ubavg(i,j,nl) enddo enddo c Nudging u c --------- margin = 0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu_w(j) do i=max(1-margin,ifu_w(j,l)),min(ii+margin,ilu_w(j,l)) ia=i-1 if (iu(ia,j).ne.1) ia=i ib=i+1 if (iu(ib,j).ne.1) ib=i ja=j-1 if (iu(i,ja).ne.1) ja=j jb=j+1 if (iu(i,jb).ne.1) jb=j ubavg(i,j,nl) = ( util1(i,j)*util5(i,j)+ . (hcu(i,j)-util1(i,j))*(util1(ia,j)*util5(ia,j)+ . util1(ib,j)*util5(ib,j)+ . util1(i,ja)*util5(i,ja)+ . util1(i,jb)*util5(i,jb))/ . (epsu(i,j)+util1(ia,j)+util1(ib,j)+ . util1(i,ja)+util1(i,jb)) )*invhcu(i,j) enddo enddo enddo c margin = 1 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) util1(i,j)=max(0., . min(util3(i,j-1),util3(i,j),hcv(i,j))) enddo i = max(1-margin,ifv_w(j,l)-1) util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcv(i,j))) if (iv_nw(i,j).eq.1) util6(i,j)=vbavg(i,j,nl) i = min(ii+margin,ilv_w(j,l)+1) util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcv(i,j))) if (iv_nw(i,j).eq.1) util6(i,j)=vbavg(i,j,nl) enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do i=1-margin,ii+margin !!Alex bug corrected from jj to ii do l=1,jsv_w(i) j=max(1-margin,jfv_w(i,l)-1) util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcv(i,j))) if (iv_nw(i,j).eq.1) util6(i,j)=vbavg(i,j,nl) j=min(jj+margin,jlv_w(i,l)+1) util1(i,j)=max(0.,min(util3(i,j-1),util3(i,j),hcv(i,j))) if (iv_nw(i,j).eq.1) util6(i,j)=vbavg(i,j,nl) enddo enddo c Nudging v c --------- margin = 0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv_w(j) do i=max(1-margin,ifv_w(j,l)),min(ii+margin,ilv_w(j,l)) ia=i-1 if (iv(ia,j).ne.1) ia=i ib=i+1 if (iv(ib,j).ne.1) ib=i ja=j-1 if (iv(i,ja).ne.1) ja=j jb=j+1 if (iv(i,jb).ne.1) jb=j vbavg(i,j,nl) = ( util1(i,j)*util6(i,j)+ . (hcv(i,j)-util1(i,j))*(util1(ia,j)*util6(ia,j)+ . util1(ib,j)*util6(ib,j)+ . util1(i,ja)*util6(i,ja)+ . util1(i,jb)*util6(i,jb))/ . (epsv(i,j)+util1(ia,j)+util1(ib,j)+ . util1(i,ja)+util1(i,jb)) )*invhcv(i,j) enddo enddo enddo c wblpf = wbaro c 840 continue ! lll=1,lstep1 c do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) uflxba(i,j) = dlt*uflxba(i,j)/delt1 enddo enddo enddo do j=1-margin,jj+margin do l=1,isu(j) i=max(1-nbdy,ifu(j,l)-1) uflxba(i,j) = dlt*uflxba(i,j)/delt1 i=min(ii+nbdy,ilu(j,l)+1) uflxba(i,j) = dlt*uflxba(i,j)/delt1 enddo enddo do j=1-margin,jj+margin do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vflxba(i,j) = dlt*vflxba(i,j)/delt1 enddo enddo enddo do i=1-margin,ii+margin do l=1,jsv(i) j=max(1-nbdy,jfv(i,l)-1) vflxba(i,j) = dlt*vflxba(i,j)/delt1 j=min(jj+nbdy,jlv(i,l)+1) vflxba(i,j) = dlt*vflxba(i,j)/delt1 enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) ubavg(i,j,n) = ubavg(i,j,nl) enddo enddo do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vbavg(i,j,n) = vbavg(i,j,nl) enddo enddo do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pbavg(i,j,n) = pbavg(i,j,nl) enddo enddo enddo c new 1 + eta c ----------- call xctilr(pbavg(1-nbdy,1-nbdy,n),1, 1, nbdy,nbdy, halo_ps) !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) oneta(i,j,n) = onetai(i,j) + pbavg(i,j,n) /pbot(i,j) enddo enddo enddo c Asselin filter c -------------- do j=1-margin,jj+margin c do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) ubavg(i,j,m) = ubavo(i,j) +0.5*ra2fac*ubavg(i,j,n) enddo enddo do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) vbavg(i,j,m) = vbavo(i,j) +0.5*ra2fac*vbavg(i,j,n) enddo enddo do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pbavg(i,j,m) = pbavo(i,j) +0.5*ra2fac*pbavg(i,j,n) enddo enddo c enddo c do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) onetao(i,j,m) = onetai(i,j) + pbavg(i,j,m) /pbot(i,j) enddo enddo enddo call xctilr(oneta( 1-nbdy,1-nbdy,n),1, 1, 2,2, halo_ps) call xctilr(onetao( 1-nbdy,1-nbdy,m),1, 1, 2,2, halo_ps) c !!Alex if (btrmas) then if (btrmas .and. .not.btrmas) then ccc if (btrmas) then !Not Validated cf. Remy Baraille c Final integration of dp c ----------------------- c call xctilr(uflxba( 1-nbdy,1-nbdy),1, 1, 2,2, halo_uv) call xctilr(vflxba( 1-nbdy,1-nbdy),1, 1, 2,2, halo_vv) call xctilr(onetacnt(1-nbdy,1-nbdy),1,1, 2,2, halo_ps) c margin = 1 util1(:,:) = 0. util2(:,:) = 0. do k = 1, kk !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) dp(i,j,k,n)=onetacnt(i,j)*dp(i,j,k,n) p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) enddo enddo enddo !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) util1(i,j) = util1(i,j) + uflx(i,j,k) enddo i=ifu(j,l)-1 if (i.ge.1-margin ) then if (iuopn(i,j).ne.0) util1(i,j)=util1(i,j) + uflx(i,j,k) endif i=ilu(j,l)+1 if (i.le.ii+margin ) then if (iuopn(i,j).ne.0) util1(i,j)=util1(i,j) + uflx(i,j,k) endif enddo enddo !$OMP END PARALLEL DO c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) util2(i,j) = util2(i,j) + vflx(i,j,k) enddo enddo enddo !$OMP END PARALLEL DO c do i=1-margin,ii+margin do l=1,jsv(i) j=jfv(i,l)-1 if (j.ge.1-margin ) then if (ivopn(i,j).ne.0) util2(i,j)=util2(i,j) + vflx(i,j,k) endif j=jlv(i,l)+1 if (j.le.jj+margin ) then if (ivopn(i,j).ne.0) util2(i,j)=util2(i,j) + vflx(i,j,k) endif enddo enddo enddo c util3(:,:) = 0. util4(:,:) = 0. do k=1,kk c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin c do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) if (uflxba(i,j)-util1(i,j).ge.0.) then if (dp(i,j,k,n).gt.0.) then util3(i,j) = util3(i,j) + dp(i-1,j,k,n) . /p(i-1,j,kk+1) endif else if (dp(i-1,j,k,n).gt.0.) then util3(i,j) = util3(i,j) + dp(i,j,k,n)/p(i,j,kk+1) endif endif enddo enddo c do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) if (vflxba(i,j)-util2(i,j).ge.0.) then if (dp(i,j,k,n).gt.0.) then util4(i,j) = util4(i,j) + dp(i,j-1,k,n) . /p(i,j-1,kk+1) endif else if (dp(i,j-1,k,n).gt.0.) then util4(i,j) = util4(i,j) + dp(i,j,k,n)/p(i,j,kk+1) endif endif enddo enddo enddo !$OMP END PARALLEL DO c enddo c do k=1,kk c !$OMP PARALLEL DO PRIVATE(j,l,i,z1) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin c do l=1,isu(j) do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) if (uflxba(i,j)-util1(i,j).ge.0.) then if (dp(i,j,k,n).gt.0.) then z1=dp(i-1,j,k,n)/p(i-1,j,kk+1)/util3(i,j) else z1=0. endif else if (dp(i-1,j,k,n).gt.0.) then z1=dp(i ,j,k,n)/p(i ,j,kk+1)/util3(i,j) else z1=0. endif endif uflux(i,j)=(uflxba(i,j)-util1(i,j))*z1 uflx(i,j,k)=uflx(i,j,k)+uflux(i,j) enddo i=ifu(j,l)-1 if (i.ge.1-margin ) then if (iuopn(i,j).ne.0) then z1=dp(i,j,k,n)/p(i,j,kk+1) uflux(i,j)=(uflxba(i,j)-util1(i,j))*z1 uflx(i,j,k)=uflx(i,j,k)+uflux(i,j) endif endif i=ilu(j,l)+1 if (i.le.ii+margin ) then if (iuopn(i,j).ne.0) then z1=dp(i-1,j,k,n)/p(i-1,j,kk+1) uflux(i,j)=(uflxba(i,j)-util1(i,j))*z1 uflx(i,j,k)=uflx(i,j,k)+uflux(i,j) endif endif enddo c do l=1,isv(j) do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) if (vflxba(i,j)-util2(i,j).ge.0.) then if (dp(i,j,k,n).gt.0.) then z1=dp(i,j-1,k,n)/p(i,j-1,kk+1)/util4(i,j) else z1=0. endif else if (dp(i,j-1,k,n).gt.0.) then z1=dp(i,j ,k,n)/p(i,j ,kk+1)/util4(i,j) else z1=0. endif endif vflux(i,j)=(vflxba(i,j)-util2(i,j))*z1 vflx(i,j,k)=vflx(i,j,k)+vflux(i,j) enddo enddo enddo !$OMP END PARALLEL DO c do i=1-margin,ii+margin do l=1,jsv(i) j=jfv(i,l)-1 if (j.ge.1-margin ) then if (ivopn(i,j).ne.0) then z1=dp(i,j,k,n)/p(i,j,kk+1) vflux(i,j)=(vflxba(i,j)-util2(i,j))*z1 vflx(i,j,k)=vflx(i,j,k)+vflux(i,j) endif endif j=jlv(i,l)+1 if (j.le.jj+margin ) then if (ivopn(i,j).ne.0) then z1=dp(i,j-1,k,n)/p(i,j-1,kk+1) vflux(i,j)=(vflxba(i,j)-util2(i,j))*z1 vflx(i,j,k)=vflx(i,j,k)+vflux(i,j) endif endif enddo enddo c !$OMP PARALLEL DO PRIVATE(j,l,i,dpmin) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) dp(i,j,k,n)=dp(i,j,k,n)- & ((uflux(i+1,j)-uflux(i,j))+ & (vflux(i,j+1)-vflux(i,j)))*delt1*scp2i(i,j) ccc p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) p(i,j,k+1)=p(i,j,k)+max(0.,dp(i,j,k,n)) enddo enddo enddo c enddo c !$OMP PARALLEL DO PRIVATE(j,l,k,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isp(j) do k=1,kk do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) if (p(i,j,kk+1) .gt. 0.) then ccc dp(i,j,k,n)=dp(i,j,k,n)*pbot(i,j)/p(i,j,kk+1) dp(i,j,k,n)=max(0.,dp(i,j,k,n))*pbot(i,j)/p(i,j,kk+1) else dp(i,j,k,n)=0. dp(i,j,1,n)=pbot(i,j) endif p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) if (isopyc .and. k.eq.1) then dpmixl(i,j,n)=dp(i,j,k,n) endif enddo enddo enddo enddo !$OMP END PARALLEL DO c call xctilr(dp(1-nbdy,1-nbdy,1,n),1,kk, 2,2, halo_ps) c endif c end subroutine barotp_dry_and_nodry