#if defined(ROW_LAND) #define SEA_P .true. #define SEA_U .true. #define SEA_V .true. #define SEA_Q .true. #elif defined(ROW_ALLSEA) #define SEA_P allip(j).or.ip(i,j).ne.0 #define SEA_U alliu(j).or.iu(i,j).ne.0 #define SEA_V alliv(j).or.iv(i,j).ne.0 #define SEA_Q alliq(j).or.iq(i,j).ne.0 #else #define SEA_P ip(i,j).ne.0 #define SEA_U iu(i,j).ne.0 #define SEA_V iv(i,j).ne.0 #define SEA_Q iq(i,j).ne.0 #endif module mod_momtum use mod_xc ! HYCOM communication interface implicit none c c --- module for momtum and related routines c private !! default is private public :: momtum_hs, momtum, momtum4 c #if defined(RELO) real, save, allocatable, dimension(:,:) :: #else real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: #endif & stress,stresx,stresy,dpmx,thkbop,oneta_mtm contains subroutine momtum_hs(m,n) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_pipe ! HYCOM debugging interface use mod_tides ! HYCOM tides #if defined(STOKES) use mod_stokes ! HYCOM Stokes drift #endif implicit none c integer m,n c c --- ----------------------------------------- c --- hydrostatic equation (and surface stress) c --- ----------------------------------------- c logical, parameter :: lpipe_momtum=.false. !usually .false. !!Alex real, parameter :: dragw_rho=0.01072d0*1026.0d0 !ice-ocean drag from CICE (2x default) real, parameter :: dragw_rho=0.00536d0*1026.0d0 !ice-ocean dragfrom CICE (default) real, parameter :: pairc=1013.0d0*100.0d0 !air pressure (Pa) real, parameter :: rgas=287.1d0 !gas constant (j/kg/k) real, parameter :: tzero=273.16d0 !celsius to kelvin offset c character text*12 c real dpdn,dpup,q,samo,simo,uimo,vimo,dpsur,psur,usur,vsur, & airt,vpmx,wndx,wndy,wind,cdw,pair,rair, & sumdp,sumth,thksur integer i,j,k,kl,l,margin,mbdy,hlstep c &,iffstep c data iffstep/0/ c save iffstep c * real*8 wtime * external wtime * real*8 wtime1(10),wtime2(20,kdm),wtimes c include 'stmt_fns.h' c #if defined(RELO) if (.not.allocated(stress)) then allocate( & stress(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & stresx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & stresy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dpmx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & thkbop(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & oneta_mtm(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)) call mem_stat_add( 6*(idm+2*nbdy)*(jdm+2*nbdy) ) stress = r_init stresx = r_init stresy = r_init dpmx = r_init thkbop = r_init oneta_mtm = 0.0 endif #endif c mbdy = 6 c c --- dp has up to date halos from cnuity or mod_hycom call xctilr(dpmixl(1-nbdy,1-nbdy, m),1, 1, 6,6, halo_ps) call xctilr(pbavg( 1-nbdy,1-nbdy,1 ),1, 2, 6,6, halo_ps) call xctilr(saln( 1-nbdy,1-nbdy,1,m),1, kk, 6,6, halo_ps) call xctilr(temp( 1-nbdy,1-nbdy,1,m),1, kk, 6,6, halo_ps) call xctilr(th3d( 1-nbdy,1-nbdy,1,m),1, kk, 6,6, halo_ps) call xctilr(oneta( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_ps) !!Alex call xctilr(oneta( 1-nbdy,1-nbdy, m),1, 1, 6,6, halo_ps) if (windf .and. iceflg.eq.2) then kl=max(nsigma,1) call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kl, 1,1, halo_uv) call xctilr(ubavg( 1-nbdy,1-nbdy, n),1, 1, 1,1, halo_uv) call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kl, 1,1, halo_vv) call xctilr(vbavg( 1-nbdy,1-nbdy, n),1, 1, 1,1, halo_vv) endif c c --- tidal forcing c if(tidflg.eq.2 .or. tidflg.eq.3) then hlstep=0 call tides_force(hlstep) endif c c --- hydrostatic equation (and surface stress) c * wtime1( 1) = wtime() c c --- rhs: th3d.m, temp.m, saln.m, p, pbavg.m c --- lhs: thstar, p, oneta, montg c margin = mbdy c !$OMP PARALLEL DO PRIVATE(j,k,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_P) then sumdp = 0.0 sumth = 0.0 do k=1,kk if (kapref.ne.0) then !thermobaric c c --- sigma-star is virtual potential density, as defined in c --- Sun et.al. (1999), 'Inclusion of thermobaricity in c --- isopycnic-coordinate ocean models', JPO 29 pp 2719-2729. c c --- use upper interface pressure in converting sigma to sigma-star. c --- to avoid density variations in layers intersected by bottom c if (kapref.gt.0) then thstar(i,j,k,1)=th3d(i,j,k,m) & +kappaf(temp(i,j,k,m), & saln(i,j,k,m), & th3d(i,j,k,m)+thbase, & 0.5*(p(i,j,k)+p(i,j,k+1)), !!!Alex & kapref) else thstar(i,j,k,1)=th3d(i,j,k,m) & +kappaf(temp(i,j,k,m), & saln(i,j,k,m), & th3d(i,j,k,m)+thbase, & 0.5*(p(i,j,k)+p(i,j,k+1)), !!!Alex & 2) thstar(i,j,k,2)=th3d(i,j,k,m) & +kappaf(temp(i,j,k,m), & saln(i,j,k,m), & th3d(i,j,k,m)+thbase, & 0.5*(p(i,j,k)+p(i,j,k+1)), !!!Alex & kapi(i,j)) #if defined(KEY_SIGMA_TRUE) thstar(i,j,k,2) = qthref - thbase - & qthref/(1.+(thstar(i,j,k,2)+thbase)*thref) #endif endif !kapref else !non-thermobaric thstar(i,j,k,1)=th3d(i,j,k,m) #if defined(KEY_SIGMA_TRUE) thstar(i,j,k,1) = qthref - thbase - & qthref/(1.+(thstar(i,j,k,1)+thbase)*thref) #endif endif !thermobaric:else c p(i,j,k+1)=p(i,j,k)+dp(i,j,k,m) c if (sshflg.ne.0) then sumth = sumth + dp(i,j,k,m)*th3d(i,j,k,m) sumdp = sumdp + dp(i,j,k,m) endif !sshflg enddo !k c if (sshflg.ne.0) then sumth = sumth / max( sumdp, onemm ) !vertical mean of th3d sumdp = sumdp*qonem * g !depth(m) * g steric(i,j) = sshgmn(i,j) + & (sshgmn(i,j) + sumdp) * & (thmean(i,j) - sumth) / & (1000.0+thbase+sumth) endif !sshflg c c --- store (1+eta) (= p_total/p_prime) in -oneta- ! oneta(i,j)=1.+pbavg(i,j,m)/p(i,j,kk+1) ! onetai(i,j)=oneta(i,j)-pbavg(i,j,m)/pbot(i,j) if (btrmas) then oneta_mtm(i,j)= oneta(i,j,m) else oneta_mtm(i,j)=1.0 !disable oneta_mtm endif c c --- m_prime in lowest layer: montg(i,j,kk,1)=psikk(i,j,1)+ & ( p(i,j,kk+1)*(thkk(i,j,1)-thstar(i,j,kk,1)) & -pbavg(i,j,m)*thstar(i,j,kk,1) )*thref**2 if (kapref.eq.-1) then montg(i,j,kk,2)=psikk(i,j,2)+ & ( p(i,j,kk+1)*(thkk(i,j,2)-thstar(i,j,kk,2)) & -pbavg(i,j,m)*thstar(i,j,kk,2) )*thref**2 endif !kapref.eq.-1 c c --- m_prime in remaining layers: do k=kk-1,1,-1 montg(i,j,k,1)=montg(i,j,k+1,1)+p(i,j,k+1) & *oneta_mtm(i,j) & *(thstar(i,j,k+1,1)-thstar(i,j,k,1))*thref**2 if (kapref.eq.-1) then montg(i,j,k,2)=montg(i,j,k+1,2)+p(i,j,k+1) & *oneta_mtm(i,j) & *(thstar(i,j,k+1,2)-thstar(i,j,k,2))*thref**2 endif !kapref.eq.-1 enddo !k c c --- srfhgt (used diagnostically, in mxmyaij and for tidal SAL). if (kapref.ne.-1) then montg1(i,j) = montg(i,j,1,1) else montg1(i,j) = skap(i,j) *montg(i,j,1,1) + & (1.0-skap(i,j))*montg(i,j,1,2) endif !kapref srfhgt(i,j) = montg1(i,j) + thref*pbavg(i,j,m) c cdiag if (sshflg.ne.0) then cdiag if (itest.gt.0 .and. jtest.gt.0) then cdiag write (lp,'(i9,2i5,3x,a,2f12.6,f12.2)') cdiag. nstep,itest+i0,jtest+j0, cdiag. 'sssh =', cdiag. steric(i,j),sshgmn(i,j),sumdp cdiag write (lp,'(i9,2i5,3x,a,3f12.6)') cdiag. nstep,itest+i0,jtest+j0, cdiag. 'thmn =', cdiag. sumth,thmean(i,j),1000.0+thbase+sumth cdiag write (lp,'(i9,2i5,3x,a,3f12.6)') cdiag. nstep,itest+i0,jtest+j0, cdiag. 'ssh =', cdiag. srfhgt(i,j),steric(i,j),srfhgt(i,j)-steric(i,j) cdiag endif !test cdiag endif !sshflg endif !ip enddo !i enddo !j !$OMP END PARALLEL DO c if (lpipe .and. lpipe_momtum .and. mslprf) then c --- compare two model runs. util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l0) write (text,'(a9,i3)') 'mslprs l=',l0 call pipe_compare(util4, ip,text) util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l1) write (text,'(a9,i3)') 'mslprs l=',l1 call pipe_compare(util4, ip,text) c util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l0)-mslprs(1:ii,0:jj-1,l0) write (text,'(a9,i3)') 'mslprY l=',l0 call pipe_compare(util4, iv,text) util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l1)-mslprs(1:ii,0:jj-1,l1) write (text,'(a9,i3)') 'mslprY l=',l1 call pipe_compare(util4, iv,text) c util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l0)-mslprs(0:ii-1,1:jj,l0) write (text,'(a9,i3)') 'mslprX l=',l0 call pipe_compare(util4, iu,text) util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l1)-mslprs(0:ii-1,1:jj,l1) write (text,'(a9,i3)') 'mslprX l=',l1 call pipe_compare(util4, iu,text) endif c c call dpudpv(dpu(1-nbdy,1-nbdy,1,m), c & dpv(1-nbdy,1-nbdy,1,m), c & p,depthu,depthv, max(0,margin-1)) cCL/RB MPI bug correction 2011-01 c call xctilr(dpu( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_us) c call xctilr(dpv( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_vs) c c --- account for temporal smoothing of mid-time dpmixl. calculate the vertical c --- excursions of the coordinates immediately above and below the mixed c --- layer base, then vertically interpolate this motion to dpmixl(i,j,m) c if(hybrid .and. mxlkta) then c c --- rhs: dp, dpmixl.m c --- lhs: util1, util2 c margin = mbdy c !$OMP PARALLEL DO PRIVATE(j,i,k,dpup,dpdn,q) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_P) then util1(i,j)=0.0 util2(i,j)=0.0 do k=1,kk util1(i,j)=util2(i,j) util2(i,j)=util2(i,j)+dpo(i,j,k,m) if (util2(i,j).ge.dpmixl(i,j,m).and. & util1(i,j).lt.dpmixl(i,j,m) ) then dpup=p(i,j,k )-util1(i,j) dpdn=p(i,j,k+1)-util2(i,j) q=(util2(i,j)-dpmixl(i,j,m))/max(onemm,dpo(i,j,k,m)) dpmixl(i,j,m)=dpmixl(i,j,m)+(dpdn+q*(dpup-dpdn)) endif enddo !k endif !ip enddo !l enddo !i !$OMP END PARALLEL DO endif c c --- -------------- c --- surface stress c --- -------------- c if (windf) then margin =0 c !$OMP PARALLEL DO PRIVATE(j,i,k, !$OMP& dpsur,psur,usur,vsur,thksur, !$OMP& airt,vpmx,wndx,wndy,wind,cdw,rair, !$OMP& uimo,vimo,simo,samo) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_P) then if (wndflg.eq.4 .or. wndflg.eq.5 .or. & (iceflg.eq.2 .and. si_c(i,j).gt.0.0)) then c --- average currents over top thkcdw meters thksur = onem*min( thkcdw, depths(i,j) ) usur = 0.0 vsur = 0.0 psur = 0.0 do k= 1,kk dpsur = min( dpo(i,j,k,n), max( 0.0, thksur-psur ) ) #if defined(STOKES) usur = usur + dpsur*(u(i,j,k,n)+u(i+1,j,k,n)+ & usd(i,j,k)+usd(i+1,j,k)) vsur = vsur + dpsur*(v(i,j,k,n)+v(i,j+1,k,n)+ & vsd(i,j,k)+vsd(i,j+1,k)) #else usur = usur + dpsur*(u(i,j,k,n)+u(i+1,j,k,n)) vsur = vsur + dpsur*(v(i,j,k,n)+v(i,j+1,k,n)) #endif psur = psur + dpsur if (psur.ge.thksur) then exit endif enddo !k usur = 0.5*( usur/psur + ubavg(i, j,n) + & ubavg(i+1,j,n) ) vsur = 0.5*( vsur/psur + vbavg(i,j, n) + & vbavg(i,j+1,n) ) endif !usur,vsur c if (wndflg.eq.2 .or. wndflg.eq.3) then ! tau on p grid if (natm.eq.2) then surtx(i,j) = taux(i,j,l0)*w0+taux(i,j,l1)*w1 surty(i,j) = tauy(i,j,l0)*w0+tauy(i,j,l1)*w1 else surtx(i,j) = taux(i,j,l0)*w0+taux(i,j,l1)*w1 & +taux(i,j,l2)*w2+taux(i,j,l3)*w3 surty(i,j) = tauy(i,j,l0)*w0+tauy(i,j,l1)*w1 & +tauy(i,j,l2)*w2+tauy(i,j,l3)*w3 endif !natm elseif (wndflg.eq.1) then ! tau on u&v grids - NOT RECOMMEDED if (natm.eq.2) then surtx(i,j) = ( (taux(i,j,l0)+taux(i+1,j,l0))*w0 & +(taux(i,j,l1)+taux(i+1,j,l1))*w1)*0.5 surty(i,j) = ( (tauy(i,j,l0)+tauy(i,j+1,l0))*w0 & +(tauy(i,j,l1)+tauy(i,j+1,l1))*w1)*0.5 else surtx(i,j) = ( (taux(i,j,l0)+taux(i+1,j,l0))*w0 & +(taux(i,j,l1)+taux(i+1,j,l1))*w1 & +(taux(i,j,l2)+taux(i+1,j,l2))*w2 & +(taux(i,j,l3)+taux(i+1,j,l3))*w3)*0.5 surty(i,j) = ( (tauy(i,j,l0)+tauy(i,j+1,l0))*w0 & +(tauy(i,j,l1)+tauy(i,j+1,l1))*w1 & +(tauy(i,j,l2)+tauy(i,j+1,l2))*w2 & +(tauy(i,j,l3)+tauy(i,j+1,l3))*w3)*0.5 endif !natm else !wndflg.eq.4,5,6 c --- calculate stress from 10m winds using cd_coare or cd_core2 if (natm.eq.2) then airt = airtmp(i,j,l0)*w0+airtmp(i,j,l1)*w1 vpmx = vapmix(i,j,l0)*w0+vapmix(i,j,l1)*w1 wndx = taux(i,j,l0)*w0+ taux(i,j,l1)*w1 wndy = tauy(i,j,l0)*w0+ tauy(i,j,l1)*w1 else airt = airtmp(i,j,l0)*w0+airtmp(i,j,l1)*w1 & +airtmp(i,j,l2)*w2+airtmp(i,j,l3)*w3 vpmx = vapmix(i,j,l0)*w0+vapmix(i,j,l1)*w1 & +vapmix(i,j,l2)*w2+vapmix(i,j,l3)*w3 wndx = taux(i,j,l0)*w0+ taux(i,j,l1)*w1 & +taux(i,j,l2)*w2+ taux(i,j,l3)*w3 wndy = tauy(i,j,l0)*w0+ tauy(i,j,l1)*w1 & +tauy(i,j,l2)*w2+ tauy(i,j,l3)*w3 endif !natm wind = sqrt( wndx**2 + wndy**2 ) if (mslprf .or. flxflg.eq.6) then if (natm.eq.2) then pair = mslprs(i,j,l0)*w0+mslprs(i,j,l1)*w1 & +prsbas else pair = mslprs(i,j,l0)*w0+mslprs(i,j,l1)*w1 & +mslprs(i,j,l2)*w2+mslprs(i,j,l3)*w3 & +prsbas endif !natm else pair = pairc endif if (flxflg.eq.6) then !use virtual temperature rair = pair / (rgas * (tzero+airt) * (1.0+0.608*vpmx)) else rair = pair / (rgas * (tzero+airt)) endif c if (wndflg.eq.4 .and. flxflg.eq.6) then c --- use wind-current in place of wind for everything samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) cdw = 1.0e-3*cd_coarep(samo,vpmx,airt,pair, & temp(i,j,1,n)) surtx( i,j) = rair*cdw*samo*(wndx-usur) surty( i,j) = rair*cdw*samo*(wndy-vsur) wndocn(i,j) = samo !save for thermf elseif (wndflg.eq.4) then cdw = 1.0e-3*cd_coare(wind,vpmx,airt, & temp(i,j,1,n)) c --- use wind-current magnitude and direction for stress samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) surtx(i,j) = rair*cdw*samo*(wndx-usur) surty(i,j) = rair*cdw*samo*(wndy-vsur) else ! wndflg.eq.5 cdw = 1.0e-3*cd_core2(wind,vpmx,airt, & temp(i,j,1,n)) c --- use wind-current magnitude and direction for stress samo = sqrt( (wndx-usur)**2 + (wndy-vsur)**2 ) surtx(i,j) = rair*cdw*samo*(wndx-usur) surty(i,j) = rair*cdw*samo*(wndy-vsur) endif endif !wndflg c surtx(i,j) = surtx(i,j) + oftaux(i,j) surty(i,j) = surty(i,j) + oftauy(i,j) c if (ishlf(i,j).eq.0) then !under an ice shelf surtx(i,j) = 0.0 surty(i,j) = 0.0 elseif (iceflg.eq.2 .and. si_c(i,j).gt.0.0) then c --- allow for ice-ocean stress uimo = si_u(i,j) - usur vimo = si_v(i,j) - vsur simo = sqrt( uimo**2 + vimo**2 ) surtx(i,j) = (1.0-si_c(i,j))*surtx(i,j) + & si_c(i,j) *dragw_rho*simo*uimo surty(i,j) = (1.0-si_c(i,j))*surty(i,j) + & si_c(i,j) *dragw_rho*simo*vimo endif !ice-ocean stress endif !ip enddo !i enddo !j !$OMP END PARALLEL DO c call xctilr(surtx,1,1, 6,6, halo_pv) call xctilr(surty,1,1, 6,6, halo_pv) endif !windf c return end subroutine momtum_hs c real function cd_coare(wind,vpmx,airt,sst) implicit none c real wind,vpmx,airt,sst c c --- Wind stress drag coefficient * 10^3 from an approximation c --- to the COARE 3.0 bulk algorithm (Fairall et al. 2003). c c --- wind = wind speed (m/s) c --- vpmx = water vapor mixing ratio (kg/kg) c --- airt = air temperature (C) c --- sst = sea temperature (C) c c --- Ta-Ts c --- ============== c --- STABLE: 7 to 0.75 degC c --- NEUTRAL: 0.75 to -0.75 degC c --- UNSTABLE: -0.75 to -8 degC c c --- Va c --- ============== c --- Low: 1 to 5 m/s c --- High: 5 to 34 m/s c c --- vamax of 34 m/s from Sraj et al, 2013 (MWR-D-12-00228.1). c real tamts,q,qva,va c real, parameter :: vamin= 1.0, vamax=34.0 real, parameter :: tdmin= -8.0, tdmax= 7.0 real, parameter :: tzero=273.16 c real, parameter :: & as0_00=-0.06695, as0_10= 0.09966, as0_20=-0.02477, & as0_01= 0.3133, as0_11=-2.116, as0_21= 0.2726, & as0_02=-0.001473, as0_12= 4.626, as0_22=-0.5558, & as0_03=-0.004056, as0_13=-2.680, as0_23= 0.3139 real, parameter :: & as5_00= 0.55815, as5_10=-0.005593, as5_20= 0.0006024, & as5_01= 0.08174, as5_11= 0.2096, as5_21=-0.02629, & as5_02=-0.0004472, as5_12=-8.634, as5_22= 0.2121, & as5_03= 2.666e-6, as5_13= 18.63, as5_23= 0.7755 real, parameter :: & au0_00= 1.891, au0_10=-0.006304, au0_20= 0.0004406, & au0_01=-0.7182, au0_11=-0.3028, au0_21=-0.01769, & au0_02= 0.1975, au0_12= 0.3120, au0_22= 0.01303, & au0_03=-0.01790, au0_13=-0.1210, au0_23=-0.003394 real, parameter :: & au5_00= 0.6497, au5_10= 0.003827, au5_20=-4.83e-5, & au5_01= 0.06993, au5_11=-0.2756, au5_21= 0.007710, & au5_02= 3.541e-5, au5_12=-1.091, au5_22=-0.2555, & au5_03=-3.428e-6, au5_13= 4.946, au5_23= 0.7654 real, parameter :: & an0_00= 1.057, an5_00= 0.6825, & an0_01=-0.06949, an5_01= 0.06945, & an0_02= 0.01271, an5_02=-0.0001029 real, parameter :: & ap0_10= as0_00 + as0_10*0.75 + as0_20*0.75**2, & ap0_11= as0_11*0.75 + as0_21*0.75**2, & ap0_12= as0_12*0.75 + as0_22*0.75**2, & ap0_13= as0_13*0.75 + as0_23*0.75**2 real, parameter :: & ap5_10= as5_00 + as5_10*0.75 + as5_20*0.75**2, & ap5_11= as5_11*0.75 + as5_21*0.75**2, & ap5_12= as5_12*0.75 + as5_22*0.75**2, & ap5_13= as5_13*0.75 + as5_23*0.75**2 real, parameter :: & am0_10= au0_00 - au0_10*0.75 + au0_20*0.75**2, & am0_11= - au0_11*0.75 + au0_21*0.75**2, & am0_12= - au0_12*0.75 + au0_22*0.75**2, & am0_13= - au0_13*0.75 + au0_23*0.75**2 real, parameter :: & am5_10= au5_00 - au5_10*0.75 + au5_20*0.75**2, & am5_11= - au5_11*0.75 + au5_21*0.75**2, & am5_12= - au5_12*0.75 + au5_22*0.75**2, & am5_13= - au5_13*0.75 + au5_23*0.75**2 c c --- saturation specific humidity (lowe, j.appl.met., 16, 100-103, 1976) real qsatur,t qsatur(t)=.622e-3*(6.107799961e+00+t*(4.436518521e-01 & +t*(1.428945805e-02+t*(2.650648471e-04 & +t*(3.031240396e-06+t*(2.034080948e-08 & +t* 6.136820929e-11)))))) c c --- correct tamts to 100% humidity tamts = airt-sst - 0.61*(airt+tzero)*(qsatur(airt)-vpmx) tamts = min( tdmax, max( tdmin, tamts ) ) va = max( vamin, min( vamax, wind ) ) qva = 1.0/va if (va.le.5.0) then if (tamts.ge. 0.75) then cd_coare = & (as0_00 + as0_01* va + as0_02* va**2 + as0_03* va**3) & + (as0_10 + as0_11*qva + as0_12*qva**2 + as0_13*qva**3) & *tamts & + (as0_20 + as0_21*qva + as0_22*qva**2 + as0_23*qva**3) & *tamts**2 elseif (tamts.le.-0.75) then cd_coare = & (au0_00 + au0_01* va + au0_02* va**2 + au0_03* va**3) & + (au0_10 + au0_11*qva + au0_12*qva**2 + au0_13*qva**3) & *tamts & + (au0_20 + au0_21*qva + au0_22*qva**2 + au0_23*qva**3) & *tamts**2 elseif (tamts.ge. -0.098) then q = (tamts+0.098)/0.848 !linear between 0.75 and -0.098 cd_coare = q* & ( ( as0_01* va + as0_02* va**2 + as0_03* va**3) & + (ap0_10 + ap0_11*qva + ap0_12*qva**2 + ap0_13*qva**3) & ) + (1.0-q)* & (an0_00 + an0_01* va + an0_02* va**2) else q = (-tamts-0.098)/0.652 !linear between -0.75 and -0.098 cd_coare = q* & ( ( au0_01* va + au0_02* va**2 + au0_03* va**3) & + (am0_10 + am0_11*qva + am0_12*qva**2 + am0_13*qva**3) & ) + (1.0-q)* & (an0_00 + an0_01* va + an0_02* va**2) endif !tamts else !va>5 if (tamts.ge. 0.75) then cd_coare = & (as5_00 + as5_01* va + as5_02* va**2 + as5_03* va**3) & + (as5_10 + as5_11*qva + as5_12*qva**2 + as5_13*qva**3) & *tamts & + (as5_20 + as5_21*qva + as5_22*qva**2 + as5_23*qva**3) & *tamts**2 elseif (tamts.le.-0.75) then cd_coare = & (au5_00 + au5_01* va + au5_02* va**2 + au5_03* va**3) & + (au5_10 + au5_11*qva + au5_12*qva**2 + au5_13*qva**3) & *tamts & + (au5_20 + au5_21*qva + au5_22*qva**2 + au5_23*qva**3) & *tamts**2 elseif (tamts.ge. -0.098) then q = (tamts+0.098)/0.848 !linear between 0.75 and -0.098 cd_coare = q* & ( ( as5_01* va + as5_02* va**2 + as5_03* va**3) & + (ap5_10 + ap5_11*qva + ap5_12*qva**2 + ap5_13*qva**3) & ) + (1.0-q)* & (an5_00 + an5_01* va + an5_02* va**2) else q = (-tamts-0.098)/0.652 !linear between -0.75 and -0.098 cd_coare = q* & ( ( au5_01* va + au5_02* va**2 + au5_03* va**3) & + (am5_10 + am5_11*qva + am5_12*qva**2 + am5_13*qva**3) & ) + (1.0-q)* & (an5_00 + an5_01* va + an5_02* va**2) endif !tamts endif !va c end function cd_coare c real function cd_coarep(samo,vpmx,airt,pair,sst) implicit none c real samo,vpmx,airt,pair,sst c c --- Wind stress drag coefficient * 10^3 from an approximation c --- to the COARE 3.0 bulk algorithm (Fairall et al. 2003). c c --- samo = wind-ocean speed (m/s) c --- vpmx = water vapor mixing ratio (kg/kg) c --- airt = air temperature (C) c --- pair = air pressure (Pa) c --- sst = sea temperature (C) c c --- Ta-Ts c --- ============== c --- STABLE: 7 to 0.75 degC c --- NEUTRAL: 0.75 to -0.75 degC c --- UNSTABLE: -0.75 to -8 degC c c --- Va c --- ============== c --- Low: 1 to 5 m/s c --- High: 5 to 34 m/s c c --- vamax of 34 m/s from Sraj et al, 2013 (MWR-D-12-00228.1). c real tamts,q,qva,va c real, parameter :: vamin= 1.0, vamax=34.0 real, parameter :: tdmin= -8.0, tdmax= 7.0 real, parameter :: tzero=273.16 c real, parameter :: & as0_00=-0.06695, as0_10= 0.09966, as0_20=-0.02477, & as0_01= 0.3133, as0_11=-2.116, as0_21= 0.2726, & as0_02=-0.001473, as0_12= 4.626, as0_22=-0.5558, & as0_03=-0.004056, as0_13=-2.680, as0_23= 0.3139 real, parameter :: & as5_00= 0.55815, as5_10=-0.005593, as5_20= 0.0006024, & as5_01= 0.08174, as5_11= 0.2096, as5_21=-0.02629, & as5_02=-0.0004472, as5_12=-8.634, as5_22= 0.2121, & as5_03= 2.666e-6, as5_13= 18.63, as5_23= 0.7755 real, parameter :: & au0_00= 1.891, au0_10=-0.006304, au0_20= 0.0004406, & au0_01=-0.7182, au0_11=-0.3028, au0_21=-0.01769, & au0_02= 0.1975, au0_12= 0.3120, au0_22= 0.01303, & au0_03=-0.01790, au0_13=-0.1210, au0_23=-0.003394 real, parameter :: & au5_00= 0.6497, au5_10= 0.003827, au5_20=-4.83e-5, & au5_01= 0.06993, au5_11=-0.2756, au5_21= 0.007710, & au5_02= 3.541e-5, au5_12=-1.091, au5_22=-0.2555, & au5_03=-3.428e-6, au5_13= 4.946, au5_23= 0.7654 real, parameter :: & an0_00= 1.057, an5_00= 0.6825, & an0_01=-0.06949, an5_01= 0.06945, & an0_02= 0.01271, an5_02=-0.0001029 real, parameter :: & ap0_10= as0_00 + as0_10*0.75 + as0_20*0.75**2, & ap0_11= as0_11*0.75 + as0_21*0.75**2, & ap0_12= as0_12*0.75 + as0_22*0.75**2, & ap0_13= as0_13*0.75 + as0_23*0.75**2 real, parameter :: & ap5_10= as5_00 + as5_10*0.75 + as5_20*0.75**2, & ap5_11= as5_11*0.75 + as5_21*0.75**2, & ap5_12= as5_12*0.75 + as5_22*0.75**2, & ap5_13= as5_13*0.75 + as5_23*0.75**2 real, parameter :: & am0_10= au0_00 - au0_10*0.75 + au0_20*0.75**2, & am0_11= - au0_11*0.75 + au0_21*0.75**2, & am0_12= - au0_12*0.75 + au0_22*0.75**2, & am0_13= - au0_13*0.75 + au0_23*0.75**2 real, parameter :: & am5_10= au5_00 - au5_10*0.75 + au5_20*0.75**2, & am5_11= - au5_11*0.75 + au5_21*0.75**2, & am5_12= - au5_12*0.75 + au5_22*0.75**2, & am5_13= - au5_13*0.75 + au5_23*0.75**2 c real satvpr,qsaturp,t,t6,p6 c c --- saturation vapor pressure (Pa), c --- from a polynominal approximation (lowe, j.appl.met., 16, 100-103, 1976) satvpr(t)= 100.0*(6.107799961e+00+t*(4.436518521e-01 & +t*(1.428945805e-02+t*(2.650648471e-04 & +t*(3.031240396e-06+t*(2.034080948e-08 & +t* 6.136820929e-11)))))) c c --- pressure dependent saturation mixing ratio (kg/kg) c --- p6 is pressure in Pa qsaturp(t6,p6)=0.622*(satvpr(t6)/(p6-satvpr(t6))) c c --- correct tamts to 100% humidity tamts = airt-sst - & 0.608*(airt+tzero)*(qsaturp(airt,pair)-vpmx) tamts = min( tdmax, max( tdmin, tamts ) ) va = max( vamin, min( vamax, samo ) ) qva = 1.0/va if (va.le.5.0) then if (tamts.ge. 0.75) then cd_coarep = & (as0_00 + as0_01* va + as0_02* va**2 + as0_03* va**3) & + (as0_10 + as0_11*qva + as0_12*qva**2 + as0_13*qva**3) & *tamts & + (as0_20 + as0_21*qva + as0_22*qva**2 + as0_23*qva**3) & *tamts**2 elseif (tamts.le.-0.75) then cd_coarep = & (au0_00 + au0_01* va + au0_02* va**2 + au0_03* va**3) & + (au0_10 + au0_11*qva + au0_12*qva**2 + au0_13*qva**3) & *tamts & + (au0_20 + au0_21*qva + au0_22*qva**2 + au0_23*qva**3) & *tamts**2 elseif (tamts.ge. -0.098) then q = (tamts+0.098)/0.848 !linear between 0.75 and -0.098 cd_coarep = q* & ( ( as0_01* va + as0_02* va**2 + as0_03* va**3) & + (ap0_10 + ap0_11*qva + ap0_12*qva**2 + ap0_13*qva**3) & ) + (1.0-q)* & (an0_00 + an0_01* va + an0_02* va**2) else q = (-tamts-0.098)/0.652 !linear between -0.75 and -0.098 cd_coarep = q* & ( ( au0_01* va + au0_02* va**2 + au0_03* va**3) & + (am0_10 + am0_11*qva + am0_12*qva**2 + am0_13*qva**3) & ) + (1.0-q)* & (an0_00 + an0_01* va + an0_02* va**2) endif !tamts else !va>5 if (tamts.ge. 0.75) then cd_coarep = & (as5_00 + as5_01* va + as5_02* va**2 + as5_03* va**3) & + (as5_10 + as5_11*qva + as5_12*qva**2 + as5_13*qva**3) & *tamts & + (as5_20 + as5_21*qva + as5_22*qva**2 + as5_23*qva**3) & *tamts**2 elseif (tamts.le.-0.75) then cd_coarep = & (au5_00 + au5_01* va + au5_02* va**2 + au5_03* va**3) & + (au5_10 + au5_11*qva + au5_12*qva**2 + au5_13*qva**3) & *tamts & + (au5_20 + au5_21*qva + au5_22*qva**2 + au5_23*qva**3) & *tamts**2 elseif (tamts.ge. -0.098) then q = (tamts+0.098)/0.848 !linear between 0.75 and -0.098 cd_coarep = q* & ( ( as5_01* va + as5_02* va**2 + as5_03* va**3) & + (ap5_10 + ap5_11*qva + ap5_12*qva**2 + ap5_13*qva**3) & ) + (1.0-q)* & (an5_00 + an5_01* va + an5_02* va**2) else q = (-tamts-0.098)/0.652 !linear between -0.75 and -0.098 cd_coarep = q* & ( ( au5_01* va + au5_02* va**2 + au5_03* va**3) & + (am5_10 + am5_11*qva + am5_12*qva**2 + am5_13*qva**3) & ) + (1.0-q)* & (an5_00 + an5_01* va + an5_02* va**2) endif !tamts endif !va c end function cd_coarep c real function cd_core2(wind,vpmx,airt,sst) implicit none c real wind,vpmx,airt,sst c c --- Wind stress drag coefficient * 10^3 from c --- the CORE v2 bulk algorithm (Large and Yeager, 2009). c c --- wind = wind speed (m/s) c --- vpmx = water vapor mixing ratio (kg/kg) c --- airt = air temperature (C) c --- sst = sea temperature (C) c c --- Added to HYCOM by Alexandra Bozec, FSU. c integer it_a real u10,v10,uw10,uw real cd_n10,cd_n10_rt,ce_n10,ch_n10,cd_rt,stab, & tv,tstar,qstar,bstar,zeta,x2,x,xx, & psi_m,psi_h,z0,rair,qrair,zi real cd10,ce10,ch10,ustar c real, parameter :: vonkar= 0.4 !Von Karmann constant real, parameter :: tzero=273.16 !celsius to kelvin offset real, parameter :: g= 9.806 !same as HYCOM's g c c --- saturation specific humidity real qsatur5,t,qra qsatur5(t,qra)= 0.98*qra*6.40380e5*exp(-5107.4/(t+tzero)) c c --- CORE v2 Large and Yeager 2009 Clim. Dyn.: The global climatology c --- of an interannually varying air-sea flux dataset. c --- The bulk formulae effectively transform the problem of specifying c --- the turbulent surface fluxes (at zi=10m) into one of describing c --- the near surface atmospheric state (wind, temperature and humidity). c * write(6,'(a,1p4g14.5)') * & ' wind,vpmx,airt,sst =',wind,vpmx,airt,sst * rair = 1.22 rair = 1.0/rair zi = 10.0 tv = (airt+tzero)*(1.0+0.608*vpmx) !in Kelvin uw = max(wind, 0.5) !0.5 m/s floor on wind (undocumented NCAR) uw10 = uw !first guess 10m wind cd_n10 = (2.7/uw10+0.142+0.0764*uw10)*1.0e-3 !L-Y eqn. 6a * write(6,'(a,1p4g14.5)') * & 'cd_n10 =',cd_n10, * & (2.7/uw10)*1.0e-3,0.142*1.0e-3,(0.0764*uw10)*1.0e-3 cd_n10_rt = sqrt(cd_n10) ce_n10 = 34.6 *cd_n10_rt*1.0e-3 !L-Y eqn. 6b stab = 0.5 + sign(0.5,airt-sst) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt*1.0e-3 !L-Y eqn. 6c cd10 = cd_n10 !first guess for exchange coeff's at z ch10 = ch_n10 ce10 = ce_n10 * write(6,'(a,1p4g14.5)') * & 'uw10,psi_m,cd_n10,cd10=',uw10,0.0,cd_n10,cd10 * do it_a= 1,2 !Monin-Obukhov iteration cd_rt = sqrt(cd10) ustar = cd_rt*uw !L-Y eqn. 7a tstar = (ch10/cd_rt)*(airt-sst) !L-Y eqn. 7b qstar = (ce10/cd_rt)* & (vpmx-qsatur5(sst,qrair)) !L-Y eqn. 7c bstar = g*(tstar/tv+qstar/(vpmx+1.0/0.608)) zeta = vonkar*bstar*zi/(ustar*ustar) !L-Y eqn. 8a zeta = sign( min(abs(zeta),10.0), zeta ) !undocumented NCAR x2 = sqrt(abs(1.0-16.0*zeta)) !L-Y eqn. 8b x2 = max(x2, 1.0) !undocumented NCAR x = sqrt(x2) if (zeta > 0.0) then psi_m = -5.0*zeta !L-Y eqn. 8c psi_h = -5.0*zeta !L-Y eqn. 8c else psi_m = log((1.0+2.0*x+x2)*(1.0+x2)/8.0) & - 2.0*(atan(x)-atan(1.0)) !L-Y eqn. 8d psi_h = 2.0*log((1.0+x2)/2.0) !L-Y eqn. 8e end if uw10 = uw/(1.0+cd_n10_rt*(log(zi/10.0)-psi_m) !L-Y eqn. 9 & /vonkar) cd_n10 = (2.7/uw10+0.142+0.0764*uw10)*1.0e-3 !L-Y eqn. 6a again cd_n10_rt = sqrt(cd_n10) ce_n10 = 34.6*cd_n10_rt*1.0e-3 !L-Y eqn. 6b again stab = 0.5 + sign(0.5,zeta) ch_n10 =(18.0*stab+32.7*(1.0-stab))*cd_n10_rt*1.0e-3 !L-Y eqn. 6c again z0 = 10.0*exp(-vonkar/cd_n10_rt) !diagnostic c xx = (log(zi/10.0)-psi_m)/vonkar cd10 = cd_n10/(1.0+cd_n10_rt*xx)**2 !L-Y 10a xx = (log(zi/10.0)-psi_h)/vonkar ch10 = ch_n10/(1.0+ch_n10*xx/cd_n10_rt) * & sqrt(cd10/cd_n10) !L-Y 10b ce10 = ce_n10/(1.0+ce_n10*xx/cd_n10_rt) * & sqrt(cd10/cd_n10) !L-Y 10c * write(6,'(a,1p4g14.5)') * & 'uw10,psi_m,cd_n10,cd10=',uw10,psi_m,cd_n10,cd10 enddo c cd_core2 = cd10*1.0e3 return end function cd_core2 c subroutine momtum(m,n) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_pipe ! HYCOM debugging interface use mod_tides ! HYCOM tides #if defined(STOKES) use mod_stokes ! Stokes Drift Velocity Module #endif use mod_nc !Netcdf output implicit none c integer m,n c c --- --------------------------------------------------------- c --- momentum equations (2nd order version) c c --- Enstrophy conserving advection scheme (Sadourney, 1975) c c --- diffusion is Laplacian and/or biharmonic, both with c --- "constant" and deformation dependent coefficients. c c --- hydrostatic equation and surface stress via momtum_hs c --- --------------------------------------------------------- c --- on entry: c --- saln(:,:,:,m) = time step t with RA time smoothing c --- saln(:,:,:,n) = time step t+1 c c --- dp(:,:,:,m) = time step t with RA time smoothing c --- dp(:,:,:,n) = time step t+1 c c --- dpo(:,:,:,n) = time step t-1 c --- dpo(:,:,:,m) = time step t c c --- dpv(:,:,:,n) = time step t-1 c --- dpv(:,:,:,m) = time step t c --- dpu(:,:,:,n) = time step t-1 c --- dpu(:,:,:,m) = time step t c c --- v(:,:,:,n) = time step t-1 c --- v(:,:,:,m) = time step t c --- u(:,:,:,n) = time step t-1 c --- u(:,:,:,m) = time step t c c --- on exit: c --- dpv(:,:,:,m) = time step t with RA time smoothing c --- dpv(:,:,:,n) = time step t+1 c --- dpu(:,:,:,m) = time step t with RA time smoothing c --- dpu(:,:,:,n) = time step t+1 c c --- v(:,:,:,m) = time step t with RA time smoothing c --- v(:,:,:,n) = time step t+1 c --- u(:,:,:,m) = time step t with RA time smoothing c --- u(:,:,:,n) = time step t+1 c c --- vtotn(:,:) = time step t-1 to t+1 barotropic tendency c --- utotn(:,:) = time step t-1 to t+1 barotropic tendency c --- --------------------------------------------------------- c logical, parameter :: lpipe_momtum=.false. !usually .false. c #if defined(RELO) real, save, allocatable, dimension(:,:) :: & vis2u,vis4u,vis2v,vis4v,vort, & wgtia,wgtib,wgtja,wgtjb, & dl2u,dl2uja,dl2ujb,dl2v,dl2via,dl2vib, & pnk0,pnkp,stresl #else real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & vis2u,vis4u,vis2v,vis4v,vort, & wgtia,wgtib,wgtja,wgtjb, & dl2u,dl2uja,dl2ujb,dl2v,dl2via,dl2vib, & pnk0,pnkp,stresl common/momtumr4/ & vis2u,vis4u,vis2v,vis4v,vort, & wgtia,wgtib,wgtja,wgtjb, & dl2u,dl2uja,dl2ujb,dl2v,dl2via,dl2vib, & pnk0,pnkp,stresl save /momtumr4/ #endif c #if defined(STOKES) #if defined(RELO) real, save, allocatable, dimension(:,:) :: #else real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: #endif & usdflux,vsdflux,dz,rholocall,rholocalr, & zbot,ztop c #if defined(RELO) real, save, allocatable, dimension(:,:,:) :: #else real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk) :: #endif & dvdzu,dudzu,dvdzv,dudzv,dzponddxu,dzponddyv, & Vtopced,Vbotced,Utopced,Ubotced,grho,dzdx,dzdy, & Vstopcedu,Vstopcedv,Vsbotcedu,Vsbotcedv, & Ustopcedu,Ustopcedv,Usbotcedu,Usbotcedv c integer km1, kp1 real uatvk0,uatvkm,uatvkp,usdatv,vatuk0,vatukm,vatukp,vsdatu #endif /* STOKES */ integer, save :: ifirst c real dpia,dpib,dpja,dpjb,vis2a,vis4a,vis2b,vis4b, & scuya,scuyb,scvxa,scvxb,vmag,dall,adrlim, & dpxy,ptopl,pbotl,cutoff,qcutoff,h1,q,deform,aspy2,aspx2, & dt1inv,phi,plo,pbop,pthkbl,ubot,vbot,pstres, & dmontg,dthstr,dragu,dragv,qdpu,qdpv,dpthin, & dpun,uhm,uh0,uhp,dpvn,vhm,vh0,vhp,sum_m,sum_n integer i,ia,ib,j,ja,jb,k,ka,l,mbdy,ktop,kmid,kbot,margin real wuv1, wuv2,oneta_u,oneta_v c * real*8 wtime * external wtime * real*8 wtime1(10),wtime2(20,kdm),wtimes c character text*12 integer, save, allocatable, dimension(:,:) :: & mask c real hfharm,a,b include 'stmt_fns.h' c c --- harmonic mean divided by 2 hfharm(a,b)=a*b/(a+b) c data ifirst / 0 / c #if defined(RELO) if (.not.allocated(stress)) then allocate( & stress(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & stresx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & stresy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dpmx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & thkbop(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & oneta_mtm(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)) call mem_stat_add( 6*(idm+2*nbdy)*(jdm+2*nbdy) ) stress = r_init stresx = r_init stresy = r_init dpmx = r_init thkbop = r_init oneta_mtm = r_init endif !stress if (.not.allocated(vis2u)) then allocate( & vis2u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vis4u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vis2v(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vis4v(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vort(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & wgtia(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & wgtib(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & wgtja(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & wgtjb(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2uja(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2ujb(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2v(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2via(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dl2vib(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & pnk0(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & pnkp(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & stresl(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) call mem_stat_add( 18*(idm+2*nbdy)*(jdm+2*nbdy) ) vis2u = 0.0 !r_init vis4u = 0.0 !r_init vis2v = 0.0 !r_init vis4v = 0.0 !r_init vort = 0.0 !r_init wgtia = 0.0 !r_init wgtib = 0.0 !r_init wgtja = 0.0 !r_init wgtjb = 0.0 !r_init dl2u = 0.0 !r_init dl2uja = 0.0 !r_init dl2ujb = 0.0 !r_init dl2v = 0.0 !r_init dl2via = 0.0 !r_init dl2vib = 0.0 !r_init pnk0 = 0.0 !r_init pnkp = 0.0 !r_init stresl = 0.0 !r_init endif !vis2u #if defined(STOKES) if (.not.allocated(usdflux)) then allocate( & usdflux(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vsdflux(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & dz(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & rholocall(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & rholocalr(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & ztop(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & zbot(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) call mem_stat_add( 7*(idm+2*nbdy)*(jdm+2*nbdy) ) usdflux = 0.0 !r_init vsdflux = 0.0 !r_init dz = 0.0 !r_init rholocall = 0.0 !r_init rholocalr = 0.0 !r_init rholocalr = 0.0 !r_init ztop = 0.0 !r_init zbot = 0.0 !r_init c allocate( & dvdzu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dudzu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dvdzv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dudzv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & grho(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dzdx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dzdy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dzponddxu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dzponddyv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vtopced(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vbotced(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Utopced(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Ubotced(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vstopcedu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vstopcedv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vsbotcedu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Vsbotcedv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Ustopcedu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Ustopcedv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Usbotcedu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & Usbotcedv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk) ) call mem_stat_add( 21*(idm+2*nbdy)*(jdm+2*nbdy)*kk ) dvdzu = 0.0 !r_init dudzu = 0.0 !r_init dvdzv = 0.0 !r_init dudzv = 0.0 !r_init grho = 0.0 !r_init dzdx = 0.0 !r_init dzdy = 0.0 !r_init dzponddxu = 0.0 !r_init dzponddyv = 0.0 !r_init Vtopced = 0.0 !r_init Vbotced = 0.0 !r_init Utopced = 0.0 !r_init Ubotced = 0.0 !r_init Vstopcedu = 0.0 !r_init Vstopcedv = 0.0 !r_init Vsbotcedu = 0.0 !r_init Vsbotcedv = 0.0 !r_init Ustopcedu = 0.0 !r_init Ustopcedv = 0.0 !r_init Usbotcedu = 0.0 !r_init Usbotcedv = 0.0 !r_init endif !usdflux #endif /* STOKES */ #endif /* RELO */ c mbdy = 6 c c --- dp and dpo have up to date halos from cnuity call xctilr(dpu( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_us) call xctilr(dpv( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_vs) call xctilr(u( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_uv) call xctilr(v( 1-nbdy,1-nbdy,1,1),1,2*kk, 6,6, halo_vv) call xctilr(ubavg(1-nbdy,1-nbdy,1 ),1, 2, 6,6, halo_uv) call xctilr(vbavg(1-nbdy,1-nbdy,1 ),1, 2, 6,6, halo_vv) call xctilr(uflx( 1-nbdy,1-nbdy,1 ),1, kk, 6,6, halo_uv) call xctilr(vflx( 1-nbdy,1-nbdy,1 ),1, kk, 6,6, halo_vv) #if defined(KEY_RASSELINUV) !!Alex definition of constant wuv2 = ra2fac wuv1 = 1.-2.*ra2fac #endif #if defined(STOKES) cedric - From Ody code dz(:,:)=0. dzponddxu(:,:,:)=0. dzponddyv(:,:,:)=0. ! utotmk(:,:,:)=0.D0 utopced(:,:,:)=0. ubotced(:,:,:)=0. ! vtotmk(:,:,:)=0.D0 vtopced(:,:,:)=0. vbotced(:,:,:)=0. rholocall(:,:)=0. rholocalr(:,:)=0. c c overlay d*dz* arrays for now. c the two uses of these fields are probably equivalent. c margin = mbdy-1 c !$OMP PARALLEL DO PRIVATE(j,i,k, !$OMP& vatuk0,vatukm,vatukp, !$OMP& uatvk0,uatvkm,uatvkp) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do k= 1,kk km1= max(1,k-1) kp1= min(kk,k+1) do i=1-margin,ii+margin if (SEA_U) then vatukm = 0.25*(v(i, j, km1,m)+ & v(i-1,j, km1,m)+ & v(i, j+1,km1,m)+ & v(i-1,j+1,km1,m) ) vatuk0 = 0.25*(v(i, j, k, m)+ & v(i-1,j, k, m)+ & v(i, j+1,k, m)+ & v(i-1,j+1,k, m) ) vatukp = 0.25*(v(i, j, kp1,m)+ & v(i-1,j, kp1,m)+ & v(i, j+1,kp1,m)+ & v(i-1,j+1,kp1,m) ) dudzu(i,j,k) = & ( (u(i,j,km1,m)*max(epsil,dpu(i,j,k, m))+ & u(i,j,k, m)*max(epsil,dpu(i,j,km1,m)) ) / & (max(epsil,dpu(i,j,k, m))+ & max(epsil,dpu(i,j,km1,m)) ) & -(u(i,j,k, m)*max(epsil,dpu(i,j,kp1,m))+ & u(i,j,kp1,m)*max(epsil,dpu(i,j,k, m)) ) / & (max(epsil,dpu(i,j,k, m))+ & max(epsil,dpu(i,j,kp1,m)) ) ) & / max(epsil,dpu(i,j,k,m)) dvdzu(i,j,k) = & ( (vatukm*max(epsil,dpu(i,j,k, m))+ & vatuk0*max(epsil,dpu(i,j,km1,m)) ) / & (max(epsil,dpu(i,j,k, m))+ & max(epsil,dpu(i,j,km1,m)) ) & -(vatuk0*max(epsil,dpu(i,j,kp1,m))+ & vatukp*max(epsil,dpu(i,j,k, m)) ) / & (max(epsil,dpu(i,j,k, m))+ & max(epsil,dpu(i,j,kp1,m)) ) ) & / max(epsil,dpu(i,j,k,m)) endif !iu if (SEA_V) then uatvkm = 0.25*(u(i, j, km1,m)+ & u(i, j-1,km1,m)+ & u(i+1,j, km1,m)+ & u(i+1,j-1,km1,m) ) uatvk0 = 0.25*(u(i, j, k, m)+ & u(i, j-1,k, m)+ & u(i+1,j, k, m)+ & u(i+1,j-1,k, m) ) uatvkp = 0.25*(u(i, j, kp1,m)+ & u(i, j-1,kp1,m)+ & u(i+1,j, kp1,m)+ & u(i+1,j-1,kp1,m) ) dvdzv(i,j,k) = & ( (v(i,j,km1,m)*max(epsil,dpv(i,j,k, m))+ & v(i,j,k, m)*max(epsil,dpv(i,j,km1,m)) ) / & (max(epsil,dpv(i,j,k, m))+ & max(epsil,dpv(i,j,km1,m)) ) & -(v(i,j,k, m)*max(epsil,dpv(i,j,kp1,m))+ & v(i,j,kp1,m)*max(epsil,dpv(i,j,k, m)) ) / & (max(epsil,dpv(i,j,k, m))+ & max(epsil,dpv(i,j,kp1,m)) ) ) & / max(epsil,dpv(i,j,k,m)) dudzv(i,j,k) = & ( (uatvkm*max(epsil,dpv(i,j,k, m))+ & uatvk0*max(epsil,dpv(i,j,km1,m)) ) / & (max(epsil,dpv(i,j,k, m))+ & max(epsil,dpv(i,j,km1,m)) ) & -(uatvk0*max(epsil,dpv(i,j,kp1,m))+ & uatvkp*max(epsil,dpv(i,j,k, m)) ) / & (max(epsil,dpv(i,j,k, m))+ & max(epsil,dpv(i,j,kp1,m)) ) ) & / max(epsil,dpv(i,j,k,m)) endif !iv enddo !i enddo !k enddo !j #endif c if (ifirst.eq.0) then ifirst=1 c --- setup zero fill. margin = mbdy c do j=1-margin,jj+margin do i=1-margin,ii+margin vis2u(i,j)=0.0 vis4u(i,j)=0.0 vis2v(i,j)=0.0 vis4v(i,j)=0.0 dl2u( i,j)=0.0 dl2v( i,j)=0.0 enddo !i enddo !j endif c c --- --------------------------------------- c --- hydrostatic equation and surface stress c --- --------------------------------------- c call momtum_hs(m,n) c c +++ ++++++++++++++++++ c +++ momentum equations c +++ ++++++++++++++++++ c #if defined(STOKES) * wtime1( 4) = wtime() margin = mbdy c !$OMP PARALLEL DO PRIVATE(j,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_U) then usdflux(i,j)=0. endif !iu if (SEA_V) then vsdflux(i,j)=0. endif !iv enddo !i enddo !j !$OMP END PARALLEL DO c c Calculate the sum of the layer thicknesses c ------------------------------------------ margin = 0 do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_P) then util2(i,j) = onetai(i,j) + pbavg(i,j,n)/pbot(i,j) endif !ip enddo !i enddo !j call xctilr(util2(1-nbdy,1-nbdy),1, 1, 1,1, halo_ps) c c Calculation of grho needed for ztop/zb ot calculation do k = 1, kk do j = 1-margin,jj+margin do i=1-margin,ii+margin if (SEA_P) then #if defined(KEY_SIGMA_TRUE) grho(i,j,k) = g*(qthref+thstar(i,j,k,1)+thbase) #else grho(i,j,k) = g/(thref-(thstar(i,j,k,1)+thbase) & *thref*thref) #endif endif ! ip enddo !i enddo !j enddo !k call xctilr(grho(1-nbdy,1-nbdy,1),1, kk, 1,1, halo_ps) c c Calculation of dzdx,dzdy needed in momtum -u-v- c ----------------------------------------- c do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_U) then ztop(i,j) = -(depths(i,j)-depths(i-1,j)) zbot(i,j) = -(depths(i,j)-depths(i-1,j)) endif ! iu enddo !i enddo !j c do k = kk, 1, -1 do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_U) then ztop(i,j) = ztop(i,j) + & (dp(i, j,k,n)*util2(i, j)/grho(i, j,k)- & dp(i-1,j,k,n)*util2(i-1,j)/grho(i-1,j,k)) dzdx(i,j,k) = 0.5*(ztop(i,j)+zbot(i,j)) zbot(i,j) = ztop(i,j) endif ! iu enddo !i enddo !j enddo !k c do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_V) then ztop(i,j) = -(depths(i,j)-depths(i,j-1)) zbot(i,j) = -(depths(i,j)-depths(i,j-1)) endif ! iv enddo !i enddo !j c do k = kk, 1, -1 do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_V) then ztop(i,j) = ztop(i,j) + & (dp(i,j, k,n)*util2(i,j )/grho(i,j, k)- & dp(i,j-1,k,n)*util2(i,j-1)/grho(i,j-1,k)) dzdy(i,j,k) = 0.5*(ztop(i,j)+zbot(i,j)) zbot(i,j) = ztop(i,j) endif ! iv enddo !i enddo !j enddo !k call xctilr(grho(1-nbdy,1-nbdy,1),1, kk, nbdy,nbdy, halo_ps) call xctilr(dzdx(1-nbdy,1-nbdy,1),1, kk, nbdy,nbdy, halo_uv) call xctilr(dzdy(1-nbdy,1-nbdy,1),1, kk, nbdy,nbdy, halo_vv) #endif /* STOKES */ c c --- rhs: p, u.n+, v.n+, ubavg.n+, vbavg.n+, depthv+, pvtrop+ c --- rhs: dpmixl.m+, taux+, dpu, depthu+, dpv, tauy+ c --- lhs: util1, util2, drag, ubrhs, stresx, vbrhs, stresy c *-----disp_count = disp_count + 1 stresl(:,:) = 0.0 c if (drglim.gt.0.0) then adrlim = drglim else adrlim = 0.125 endif dt1inv = 1./delt1 c margin = mbdy - 1 c !$OMP PARALLEL DO PRIVATE(j,l,i,k, !$OMP& phi,plo,pbop,ubot,vbot,vmag,dall,pstres) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin c do i=1-margin,ii+margin if (SEA_P) then c c --- bottom drag (standard bulk formula) c --- bottom stress is applied over thickness thkbot (or the c --- total depth is this is less) c thkbop(i,j)=min( thkbot*onem, p(i,j,kk+1) ) c c --- the bottom stress term is estimated using velocity averaged c --- over the bottom boundary layer. this thickness is dpbbl for c --- the kpp boundary layer; otherwise, it is thkbop c ubot=0.0 vbot=0.0 if (mxlkpp .and. bblkpp) then pthkbl=max(dpbbl(i,j),thkbop(i,j)) !thknss of bot. b.l. else pthkbl=thkbop(i,j) !thknss of bot. b.l. endif pbop=p(i,j,kk+1)-pthkbl !top of bot. b.l. phi =max(p(i,j,1),pbop) do k=1,kk plo =phi ! max(p(i,j,k),pbop) phi =max(p(i,j,k+1),pbop) ubot=ubot + (u(i,j,k,n)+u(i+1,j,k,n))*(phi-plo) vbot=vbot + (v(i,j,k,n)+v(i,j+1,k,n))*(phi-plo) *RBc Modif RB pour friction avec cl ouvertes *RBccc ubot=ubot + (u(i,j,k,n)+u(i+1,j,k,n))*(phi-plo) *RBccc vbot=vbot + (v(i,j,k,n)+v(i,j+1,k,n))*(phi-plo) *RB if (iuopn(i,j).ge.1) then *RB ubot=ubot + (unest(i,j,k,1)+u(i+1,j,k,n))*(phi-plo) *RB elseif (iuopn(i+1,j).ge.1) then *RB ubot=ubot + (u(i,j,k,n)+unest(i+1,j,k,1))*(phi-plo) *RB else *RB ubot=ubot + (u(i,j,k,n)+u(i+1,j,k,n))*(phi-plo) *RB endif *RB if (ivopn(i,j).ge.1) then *RB vbot=vbot + (vnest(i,j,k,1)+v(i,j+1,k,n))*(phi-plo) *RB elseif (ivopn(i,j+1).ge.1) then *RB vbot=vbot + (v(i,j,k,n)+vnest(i,j+1,k,1))*(phi-plo) *RB else *RB vbot=vbot + (v(i,j,k,n)+v(i,j+1,k,n))*(phi-plo) *RB endif enddo !k ubot=ubot/min(pthkbl,p(i,j,kk+1)) & + (ubavg(i,j,n)+ubavg(i+1,j,n)) vbot=vbot/min(pthkbl,p(i,j,kk+1)) & + (vbavg(i,j,n)+vbavg(i,j+1,n)) c c --- drag = Cb * (|v| + c.bar) c --- include 1/thkbop for the fraction of layer calculation c --- and onem for conversion from 1/dp to 1/h c vmag=0.5*sqrt(ubot**2+vbot**2) dall=cbp(i,j)*(vmag+cbarp(i,j)) !no tidal drag drag(i,j)=dall*onem/thkbop(i,j) if (mxlkpp .and. bblkpp) then ustarb(i,j)=sqrt(dall*vmag) endif c c --- tidal bottom drag, drgten.1.1 is drgscl*rh c util6(i,j)=max(0.1,thkdrg)*onem !drag applied over this thknss util5(i,j)=drgten(1,1,i,j)/min(util6(i,j)*qonem,depths(i,j)) endif !ip enddo !i c c --- store r.h.s. of barotropic u/v eqn. in -ubrhs,vbrhs- c --- time-interpolate wind stress c do i=1-margin,ii+margin if (SEA_U) then ubrhs(i,j)= & (vbavg(i ,j, m)*depthv(i ,j) & +vbavg(i ,j+1,m)*depthv(i ,j+1) & +vbavg(i-1,j, m)*depthv(i-1,j) & +vbavg(i-1,j+1,m)*depthv(i-1,j+1)) & *(pvtrop(i,j)+pvtrop(i,j+1))*.125 c if (windf) then if(hybrid .and. mxlkrt) then pstres=0.5*(dpmixl(i,j,m)+dpmixl(i-1,j,m)) else pstres=dpu(i,j,1,m) endif c --- units of surtx are N/m^2 (i.e. Pa) stresx(i,j)=(surtx(i,j)+surtx(i-1,j))*0.5*g & /(pstres*thref) else ! no taux stresx(i,j)=0. endif !windf:else endif !iu c if (SEA_V) then vbrhs(i,j)= & -(ubavg(i, j ,m)*depthu(i,j ) & +ubavg(i+1,j ,m)*depthu(i+1,j ) & +ubavg(i, j-1,m)*depthu(i,j-1) & +ubavg(i+1,j-1,m)*depthu(i+1,j-1)) & *(pvtrop(i,j)+pvtrop(i+1,j))*.125 c if (windf) then if(hybrid .and. mxlkrt) then pstres=0.5*(dpmixl(i,j,m)+dpmixl(i,j-1,m)) else pstres=dpv(i,j,1,m) endif c --- units of surty are N/m^2 (i.e. Pa) stresy(i,j)=(surty(i,j)+surty(i,j-1))*0.5*g & /(pstres*thref) else ! no tauy stresy(i,j)=0. endif !windf:else endif !iv enddo !i enddo !j !$OMP END PARALLEL DO c if (lpipe .and. lpipe_momtum) then c --- compare two model runs. write (text,'(a9,i3)') 'uba.n n=',n call pipe_compare(ubavg(1-nbdy,1-nbdy,n),iu,text) write (text,'(a9,i3)') 'vba.n n=',n call pipe_compare(vbavg(1-nbdy,1-nbdy,n),iv,text) write (text,'(a9,i3)') 'drag k=',0 call pipe_compare(drag,ip,text) endif c c --- the old momeq2.f starts here c dpthin = 0.001*onemm h1 = tenm !used in lateral weighting of hor.pres.grad. cutoff = 0.5*onem qcutoff = 1.0/cutoff c * wtime1( 5) = wtime() c c --- rhs: 0.0 c --- lhs: util1, util2 c margin = mbdy - 2 c !$OMP PARALLEL DO PRIVATE(j,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin c --- spatial weighting function for pressure gradient calculation: util1(i,j)=0.0 util2(i,j)=0.0 c --- p.1 pnkp( i,j)=0.0 enddo !i enddo !j #if defined(STOKES) cOdyOdyOdy ##if defined(KEY_STOKES_GENERIC_COORDINATE_OLD) c On cree une boucle sur k pour calculer les vitesses aux interfaces c superieure et inferieure des couches. Ce calcul faisant intervenir c les vitesses de couches sucessives, le lissage temporel utilise dans c la boucle principale sur k ne modifie pas les valeurs des vitesses utilisees. c>>>>>>>>> From Babel Fish <<<<<<<<<<<<<<<<<<<< c We created a loop on k to calculate speeds to the upper and lower layers of interfaces. c This calculation involving layers successive temporal smoothing velocities used in the c main loop on k does not modify the values of speeds used. do k=1,kk c On interpole les vitesses sur les interfaces superieure et inferieure c a partir des couches superieure (k-1), actuelle (k) et inferieure (k+1) c>>>>>>>>> From Babel Fish <<<<<<<<<<<<<<<<<<<< c It interpolates the speeds on the interfaces upper and lower from the c layers top (k-1), current (k) and lower (k + 1) if (k.NE.1 .AND. k.NE.KK) then ktop=k-1 kmid=k kbot=k+1 c Sauf pour les premiere et derniere couches - on fait l'hypothese que la c pente est egale a la couche 2 pour la premiere couche et egale a l'avant derniere couche c pour la derniere couche c>>>>>>>>> From Babel Fish <<<<<<<<<<<<<<<<<<<< c Except for the first and last is layers - l 'hypothesis that the slope is equal c to the layer 2 for the first layer and equal a l' before last layer to the last c layer else c On prend les valeurs de la couche 2 if (k.EQ.1) then ktop=1 kmid=2 kbot=3 c On prend les valeurs de l'avant derniere couche else ktop=KK-2 kmid=KK-1 kbot=KK endif endif c Boucle sur frontieres points V - j'ai repris le code dans la boucle principale c (avec les commentaires) do i=1-margin,ii+margin do l=1,jsv(i) !ok c Frontiere sud du segment l considere j=jfv(i,l)-1 c cAP if (j .ge. 1-margin) then if( ivopn(i,j) .gt. 0 ) then c cAP c to be checked ! Division par somme des ponderations - Attention a division par zero Vtopced(i,j,k)=v(i,j+1,ktop,m)*dpv(i,j,ktop,m)+ & v(i,j+1,kmid,m)*dpv(i,j,kmid,m) Vtopced(i,j,k)=Vtopced(i,j,k)/(dpv(i,j,ktop,m)+ & dpv(i,j,kmid,m)+epsil) Vbotced(i,j,k)=v(i,j+1,kbot,m)*dpv(i,j,kbot,m)+ & v(i,j+1,kmid,m)*dpv(i,j,kmid,m) Vbotced(i,j,k)=Vbotced(i,j,k)/(dpv(i,j,kbot,m)+ & dpv(i,j,kmid,m)+epsil) Utopced(i,j,k)=u(i,j+1,ktop,m)*dpu(i,j,ktop,m)+ & u(i,j+1,kmid,m)*dpu(i,j,kmid,m) Utopced(i,j,k)=Utopced(i,j,k)/(dpu(i,j,ktop,m)+ & dpu(i,j,kmid,m)+epsil) Ubotced(i,j,k)=u(i,j+1,kbot,m)*dpu(i,j,kbot,m)+ & u(i,j+1,kmid,m)*dpu(i,j,kmid,m) Ubotced(i,j,k)=Ubotced(i,j,k)/(dpu(i,j,kbot,m)+ & dpu(i,j,kmid,m)+epsil) endif endif c Frontiere nord du segment l considere j=jlv(i,l)+1 if (j .le. jj+margin) then if (ivopn(i,j) .gt. 0) then c cAP c to be checked ! NB en toute rigueur, c'est 0.5*dpu mais le 0.5 est simplifie au final ! Division par somme des ponderations - Attention a division par zero Vtopced(i,j,k)=v(i,j-1,ktop,m)*dpv(i,j,ktop,m)+ & v(i,j-1,kmid,m)*dpv(i,j,kmid,m) Vtopced(i,j,k)=Vtopced(i,j,k)/(dpv(i,j,ktop,m)+ & dpv(i,j,kmid,m)+epsil) Vbotced(i,j,k)=v(i,j-1,kbot,m)*dpv(i,j,kbot,m)+ & v(i,j-1,kmid,m)*dpv(i,j,kmid,m) Vbotced(i,j,k)=Vbotced(i,j,k)/(dpv(i,j,kbot,m)+ & dpv(i,j,kmid,m)+epsil) Utopced(i,j,k)=u(i,j-1,ktop,m)*dpu(i,j,ktop,m)+ & u(i,j-1,kmid,m)*dpu(i,j,kmid,m) Utopced(i,j,k)=utopced(i,j,k)/(dpu(i,j,ktop,m)+ & dpu(i,j,kmid,m)+epsil) Ubotced(i,j,k)=u(i,j-1,kbot,m)*dpu(i,j,kbot,m)+ & u(i,j-1,kmid,m)*dpu(i,j,kmid,m) Ubotced(i,j,k)=ubotced(i,j,k)/(dpu(i,j,kbot,m)+ & dpu(i,j,kmid,m)+epsil) endif endif enddo !l enddo !i c================================================================ c DO the East and West frontiers need to be examined ?? c---------------------------------------------------------------- c c Reste des points V - domaine principale do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_V) then ! NB en toute rigueur, c'est 0.5*dpu mais le 0.5 est simplifie au final ! Division par somme des ponderations - Attention a division par zero Vtopced(i,j,k)=v(i,j,ktop,m)*dpv(i,j,ktop,m)+ & v(i,j,kmid,m)*dpv(i,j,kmid,m) Vtopced(i,j,k)=Vtopced(i,j,k)/(dpv(i,j,ktop,m)+ & dpv(i,j,kmid,m)+epsil) Vbotced(i,j,k)=v(i,j,kbot,m)*dpv(i,j,kbot,m)+ & v(i,j,kmid,m)*dpv(i,j,kmid,m) Vbotced(i,j,k)=Vbotced(i,j,k)/(dpv(i,j,kbot,m)+ & dpv(i,j,kmid,m)+epsil) Utopced(i,j,k)=u(i,j,ktop,m)*dpu(i,j,ktop,m)+ & u(i,j,kmid,m)*dpu(i,j,kmid,m) Utopced(i,j,k)=Utopced(i,j,k)/(dpu(i,j,ktop,m)+ & dpu(i,j,kmid,m)+epsil) Ubotced(i,j,k)=u(i,j,kbot,m)*dpu(i,j,kbot,m)+ & u(i,j,kmid,m)*dpu(i,j,kmid,m) Ubotced(i,j,k)=Ubotced(i,j,k)/(dpu(i,j,kbot,m)+ & dpu(i,j,kmid,m)+epsil) endif !iv enddo !i enddo !j c Fin de la boucle k enddo !k c On cree une boucle sur k pour calculer les vitesses aux interfaces c superieure et inferieure des couches. Ce calcul faisant intervenir c les vitesses de couches sucessives, le lissage temporel utilise dans c la boucle principale sur k ne modifie pas les valeurs des vitesses utilisees. c>>>>>>>>>>>>>>>>>>>>>>>>>>>> BABEL FISH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< c We created a loop on k to calculate speeds to the upper and lower c layers of interfaces. This calculation involving layers successive c temporal smoothing speeds uses ! print *,'On passe par ici' ! STOP do k=1,kk c On interpole les vitesses sur les interfaces superieure et inferieure c a partir des couches superieure (k-1), actuelle (k) et inferieure (k+1) if (k.NE.1 .AND. k.NE.KK) then ktop=k-1 kmid=k kbot=k+1 c Sauf pour les premiere et derniere couches - on fait l'hypothese que la c pente est egale a la couche 2 pour la premiere couche et egale a l'avant derniere couche c pour la derniere couche else c On prend les valeurs de la couche 2 if (k.EQ.1) then ktop=1 kmid=2 kbot=3 c On prend les valeurs de l'avant derniere couche else ktop=KK-2 kmid=KK-1 kbot=KK endif endif c !!! ON TRAVAILLE AVEC LA VITESSE STOKES DEFINIE EN TOUT PT (voir waves.f) c Reste des points V - domaine principale do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_V) then ! NB en toute rigueur, c'est 0.5*dpu mais le 0.5 est simplifie au final ! Division par somme des ponderations - Attention a division par zero Vstopcedu(i,j,k)=usd(i,j,ktop)*dpv(i,j,ktop,m)+ & usd(i,j,kmid)*dpv(i,j,kmid,m) Vstopcedu(i,j,k)=Vstopcedu(i,j,k)/(dpv(i,j,ktop,m)+ & dpv(i,j,kmid,m)+epsil) Vsbotcedu(i,j,k)=usd(i,j,kbot)*dpv(i,j,kbot,m)+ & usd(i,j,kmid)*dpv(i,j,kmid,m) Vsbotcedu(i,j,k)=Vsbotcedu(i,j,k)/(dpv(i,j,kbot,m)+ & dpv(i,j,kmid,m)+epsil) Ustopcedv(i,j,k)=vsd(i,j,ktop)*dpu(i,j,ktop,m)+ & vsd(i,j,kmid)*dpu(i,j,kmid,m) Ustopcedv(i,j,k)=Ustopcedv(i,j,k)/(dpu(i,j,ktop,m)+ & dpu(i,j,kmid,m)+epsil) Usbotcedv(i,j,k)=vsd(i,j,kbot)*dpu(i,j,kbot,m)+ & vsd(i,j,kmid)*dpu(i,j,kmid,m) Usbotcedv(i,j,k)=Usbotcedv(i,j,k)/(dpu(i,j,kbot,m)+ & dpu(i,j,kmid,m)+epsil) endif !iv enddo !i enddo !j c Fin de la boucle k enddo !k #endif c do 9 k=1,kk #if defined(STOKES) c km1= max(1,k-1) kp1= min(kk,k+1) if (k.NE.1 .AND. k.NE.KK) then ktop=k-1 kmid=k kbot=k+1 else c On prend les valeurs de la couche 2 if (k.EQ.1) then ktop=1 kmid=2 kbot=3 c On prend les valeurs de l'avant derniere couche else ktop=KK-2 kmid=KK-1 kbot=KK endif endif #endif c c --- store total (barotropic plus baroclinic) flow at old and mid time in c --- -utotn,vtotn- and -utotm,vtotm- respectively. store minimum thickness c --- values for use in pot.vort. calculation in -dpmx-. c --- store p.k and p.k+1 for time level t+1 in pnk0,pnkp c * wtime2( 1,k) = wtime() c c --- rhs: dpmx, dp.m+ c --- lhs: dpmx c margin = mbdy - 2 c do i=1-margin,ii+margin dpmx(i,1-margin)=2.*cutoff enddo !i !$OMP PARALLEL DO PRIVATE(j,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin dpmx(i,j+1)=2.*cutoff enddo !i do i=1-margin,ii+margin if (SEA_U) then dpmx(i,j+1)=max(dpmx(i,j+1),dpo(i,j,k,m)+dpo(i-1,j,k,m)) endif !iu if (SEA_P) then pnk0(i,j)=pnkp(i,j) pnkp(i,j)=pnk0(i,j)+dp(i,j,k,n) endif !ip enddo !i enddo !j c * wtime2( 2,k) = wtime() c c --- rhs: ubavg.m, ubavg.n, dp.m+, dpu c --- lhs: utotm, utotn, uflux, dpmx, pu c * margin = mbdy - 2 !!Alex c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) !ok c i=ifu(j,l)-1 if (i.ge.1-margin .and. i.le.ii+margin) then if (iuopn(i,j).ne.0) then utotm(i,j)=u(i+1,j,k,m)+ubavg(i,j,m) utotn(i,j)=u(i+1,j,k,n)+ubavg(i,j,n) uflux(i,j)=utotm(i,j)*max(dpo(i,j,k,m),cutoff) #if defined(STOKES) usdflux(i,j)=(usd(i+1,j,k)+usdbavg(i,j))* & max(dpo(i,j,k,m),cutoff) #endif endif endif i=ilu(j,l)+1 if (i.ge.1-margin .and. i.le.ii+margin) then if (iuopn(i,j).ne.0) then utotm(i,j)=u(i-1,j,k,m)+ubavg(i,j,m) utotn(i,j)=u(i-1,j,k,n)+ubavg(i,j,n) uflux(i,j)=utotm(i,j)*max(dpo(i-1,j,k,m),cutoff) #if defined(STOKES) usdflux(i,j)=(usd(i-1,j,k)+usdbavg(i,j))* & max(dpo(i-1,j,k,m),cutoff) #endif endif endif c do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) dpmx(i,j)=max(dpmx(i,j),dpo(i,j,k,m)+dpo(i-1,j,k,m)) utotm(i,j)=u(i,j,k,m)+ubavg(i,j,m) utotn(i,j)=u(i,j,k,n)+ubavg(i,j,n) uflux(i,j)=utotm(i,j)*max(dpu(i,j,k,m),cutoff) #if defined(STOKES) usdflux(i,j)=(usd(i,j,k)+usdbavg(i,j))* & max(dpu(i,j,k,m),cutoff) #endif pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,m) c enddo !i enddo !l enddo !j !$OMP END PARALLEL DO c * wtime2( 3,k) = wtime() c c --- rhs: vbavg.m, vbavg.n, dp.m+, dpv c --- lhs: vtotm, vtotn, vflux, dpmx, pv c * margin = mbdy - 2 !!Alex c do i=1-margin,ii+margin do l=1,jsv(i) !ok j=jfv(i,l)-1 if (j.ge.1-margin .and. j.le.jj+margin) then if (ivopn(i,j).ne.0) then vtotm(i,j)=v(i,j+1,k,m)+vbavg(i,j,m) vtotn(i,j)=v(i,j+1,k,n)+vbavg(i,j,n) vflux(i,j)=vtotm(i,j)*max(dpo(i,j,k,m),cutoff) #if defined(STOKES) vsdflux(i,j)=(vsd(i,j+1,k)+vsdbavg(i,j))* & max(dpo(i,j,k,m),cutoff) #endif endif endif j=jlv(i,l)+1 if (j.ge.1-margin .and. j.le.jj+margin) then if (ivopn(i,j).gt.0) then vtotm(i,j)=v(i,j-1,k,m)+vbavg(i,j,m) vtotn(i,j)=v(i,j-1,k,n)+vbavg(i,j,n) vflux(i,j)=vtotm(i,j)*max(dpo(i,j-1,k,m),cutoff) #if defined(STOKES) vsdflux(i,j)=(vsd(i,j-1,k)+vsdbavg(i,j))* & max(dpo(i,j-1,k,m),cutoff) #endif endif endif enddo !l enddo !i c * wtime2( 4,k) = wtime() !$OMP PARALLEL DO PRIVATE(j,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin if (SEA_V) then dpmx(i ,j)=max(dpmx(i ,j),dpo(i,j,k,m)+dpo(i,j-1,k,m)) dpmx(i+1,j)=max(dpmx(i+1,j),dpo(i,j,k,m)+dpo(i,j-1,k,m)) vtotm(i,j)=v(i,j,k,m)+vbavg(i,j,m) vtotn(i,j)=v(i,j,k,n)+vbavg(i,j,n) vflux(i,j)=vtotm(i,j)*max(dpv(i,j,k,m),cutoff) #if defined(STOKES) vsdflux(i,j)=(vsd(i,j,k)+vsdbavg(i,j))* & max(dpv(i,j,k,m),cutoff) #endif pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,m) endif !iv enddo !i enddo !j c c --- define auxiliary velocity fields (via,vib,uja,ujb) to implement c --- sidewall friction along near-vertical bottom slopes. wgtja,wgtjb,wgtia, c --- wgtib indicate the extent to which a sidewall is present. c * wtime2( 5,k) = wtime() c c --- rhs: pu, depthu+, utotn+, wgtja c --- lhs: wgtja, wgtjb, uja, ujb, dl2u c --- rhs: pv, depthv+, vtotn+, wgtia c --- lhs: wgtia, wgtib, via, vib, dl2v c --- rhs: vtotm, vort+, corio+, dp.m+, dpmx+, vtotn c --- lhs: vort, potvor, defor2 c * margin = mbdy - 2 !!Alex c !$OMP PARALLEL DO PRIVATE(j,ja,jb,l,i,ia,ib,aspy2,aspx2) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin c --- assume margin to properly commute pressure torque from layer to layer c> Jul. 2000 - loop reordering and logic changes for OpenMP c> Aug. 2000 - loop 117 executed only when hybrid vertical coordinate is used c> Oct. 2001 - replaced biharm with veldf[24] and visco[24] c> Sep. 2004 - kapref selects one of three thermobaric reference states c> Jun. 2006 - split out momtum_hs, for initial ssh calculation c> Nov. 2006 - inserted volume-force for tide c> Apr. 2007 - added drglim, implicit or CFL-limited explicit bottom drag c> Apr. 2007 - added dragrh, linear tidal drag based on bottom roughness c> Apr. 2007 - btrlfr: moved [uvp]bavg assignment to barotp c> Jun. 2007 - added momtum4 c> Mar 2009 - more accurate kappaf, with potential density c> Now 2009 - several bugfixes to momtum4 c> Mar 2010 - remove anti-drag on detided velocity c> Apr 2010 - put back anti-drag on detided velocity c> Apr 2011 - cbarp(i,j) in place of cbar c> Jul 2011 - salfac(i,j) in place of tidsal c> Jul 2011 - modified momtum4 based on NERSC 2.1.03_MPI momtum_quick c> Jul 2011 - modified momtum4 to use viscp and viscq c> Jul 2011 - modified momtum4 to use local arrays for y advection c> Jul 2011 - modified momtum4 to always use v 0 -v for extrapolation c> Aug 2011 - use ra2fac for wuv1 and wuv2 (so now wuv==wts) c> Aug 2011 - replaced dpold,dpoldm with dpo c> Aug 2011 - reworked Robert-Asselin time filter c> Sep 2011 - cbp(i,j) in place of cb c> Jan 2012 - added thkcdw c> July 2012 - bugfix for bottom drag when depth < thkbot c> July 2012 - thkbop is now always based on thkbot (even for bblkpp) c> Nov. 2012 - surtx,y are halo_pv (not halo_ps) c> Nov. 2012 - wndflg=4 for calculating wind stress using cd_coare c> Nov. 2012 - oftaux,oftauy for wind stress offset c> Jan. 2013 - added thkdrg=0.0 for tidal drag in barotp.f c> Jan. 2013 - replaced dragrh with drgten.1.1 c> Jan. 2013 - added gtide for tidal body forcing c> May 2013 - removed thbase from montg.kk c> July 2013 - vamax set to 34 m/s c> Sep. 2013 - optionally added Stokes drift c> Nov. 2013 - wndflg=5 for calculating wind stress using cd_core2 c> Jan. 2014 - tv is in Kelvin (cd_core2) c> Jan. 2014 - added gslpr for atmospheric pressure forcing c> Jan. 2014 - added natm c> Apr. 2014 - added pair for atmospheric pressure forcing c> Apr. 2014 - added ice shelf logic (ishlf) c> May 2014 - use land/sea masks (e.g. ip) to skip land c> Oct 2014 - added cd_coarep for flxflg=6 c> Oct 2014 - flxflg=6 uses samo in place of wind everywhere c> Oct 2014 - flxflg=4,5 uses wind-current magnitude and direction for stress c> Oct. 2042 - dragw_rho appropriate for thkcdw=3.0, 2x default c> Apr. 2015 - moved atmospheric pressure forcing to barotp c> Apr. 2015 - moved tidal body forcing to barotp