#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 ! ! --- module for momtum and related routines ! private !! default is private public :: momtum_hs, momtum, momtum4, momtum_init ! #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, & defor1, defor2, & ! deformation components uflux1,vflux1, & ! mass fluxes potvor ! potential vorticity contains subroutine momtum_init() implicit none ! ! --- ---------------------------------------------- ! --- Initialization of arrays for momentum equation ! --- ---------------------------------------------- ! #if defined(RELO) allocate( & defor1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & defor2(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & uflux1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & vflux1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & potvor(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), & 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) ) call mem_stat_add( 10*(idm+2*nbdy)*(jdm+2*nbdy) ) #endif stress = 0.0 stresx = r_init stresy = r_init dpmx = r_init thkbop = r_init ! --- all of these should be zero on land. defor1 = 0.0 defor2 = 0.0 uflux1 = 0.0 vflux1 = 0.0 potvor = 0.0 return end subroutine momtum_init 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 ! integer m,n ! ! --- ----------------------------------------- ! --- hydrostatic equation (and surface stress) ! --- ----------------------------------------- ! logical, parameter :: lpipe_momtum=.false. !usually .false. real, parameter :: dragw_rho=0.01072d0*1026.0d0 !ice-ocean drag from CICE 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 ! character text*12 ! 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 ! &,iffstep ! data iffstep/0/ ! save iffstep ! ! real*8 wtime ! external wtime ! real*8 wtime1(10),wtime2(20,kdm),wtimes ! # include "stmt_fns.h" ! ! mbdy = 6 ! ! --- 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) 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 ! ! --- tidal forcing ! if(tidflg.eq.2 .or. tidflg.eq.3) then hlstep=0 call tides_force(hlstep) endif ! ! --- hydrostatic equation (and surface stress) ! ! wtime1( 1) = wtime() ! ! --- rhs: th3d.m, temp.m, saln.m, p, pbavg.m ! --- lhs: thstar, p, oneta, montg ! margin = mbdy ! !$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 ! ! --- sigma-star is virtual potential density, as defined in ! --- Sun et.al. (1999), 'Inclusion of thermobaricity in ! --- isopycnic-coordinate ocean models', JPO 29 pp 2719-2729. ! #if defined(KAPPAF_CENTERED) ! --- use layer centered pressure in converting sigma to sigma-star. #else ! --- use upper interface pressure in converting sigma to sigma-star. ! --- to avoid density variations in layers intersected by bottom #endif ! 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, & #if defined(KAPPAF_CENTERED) 0.5*(p(i,j,k+1)+p(i,j,k)), & #else p(i,j,k), & #endif 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, & #if defined(KAPPAF_CENTERED) 0.5*(p(i,j,k+1)+p(i,j,k)), & #else p(i,j,k), & #endif 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, & #if defined(KAPPAF_CENTERED) 0.5*(p(i,j,k+1)+p(i,j,k)), & #else p(i,j,k), & #endif kapi(i,j)) endif !kapref else !non-thermobaric thstar(i,j,k,1)=th3d(i,j,k,m) endif !thermobaric:else ! p(i,j,k+1)=p(i,j,k)+dp(i,j,k,m) ! if (sshflg.eq.1) then sumth = sumth + dp(i,j,k,m)*th3d(i,j,k,m) sumdp = sumdp + dp(i,j,k,m) endif !sshflg enddo !k ! if (sshflg.eq.1) 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 ! ! --- m_prime in lowest layer: montg(i,j,kk,1)=psikk(i,j,1)+montg_c(i,j)+ & ( p(i,j,kk+1)*(thkk(i,j,1)-thstar(i,j,kk,1)) & -pbavg(i,j,m)*thstar(i,j,kk,1) )*svref**2 if (kapref.eq.-1) then montg(i,j,kk,2)=psikk(i,j,2)+montg_c(i,j)+ & ( p(i,j,kk+1)*(thkk(i,j,2)-thstar(i,j,kk,2)) & -pbavg(i,j,m)*thstar(i,j,kk,2) )*svref**2 endif !kapref.eq.-1 ! ! --- 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(i,j,m) & *(thstar(i,j,k+1,1)-thstar(i,j,k,1))*svref**2 if (kapref.eq.-1) then montg(i,j,k,2)=montg(i,j,k+1,2)+p(i,j,k+1)*oneta(i,j,m) & *(thstar(i,j,k+1,2)-thstar(i,j,k,2))*svref**2 endif !kapref.eq.-1 enddo !k ! ! --- 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) + svref*pbavg(i,j,m) ! !diag if (sshflg.eq.1) then !diag if (itest.gt.0 .and. jtest.gt.0) then !diag write (lp,'(i9,2i5,3x,a,2f12.6,f12.2)') & !diag nstep,itest+i0,jtest+j0, & !diag 'sssh =', & !diag steric(i,j),sshgmn(i,j),sumdp !diag write (lp,'(i9,2i5,3x,a,3f12.6)') & !diag nstep,itest+i0,jtest+j0, & !diag 'thmn =', & !diag sumth,thmean(i,j),1000.0+thbase+sumth !diag write (lp,'(i9,2i5,3x,a,3f12.6)') & !diag nstep,itest+i0,jtest+j0, & !diag 'ssh =', & !diag srfhgt(i,j),steric(i,j),srfhgt(i,j)-steric(i,j) !diag endif !test !diag endif !sshflg endif !ip enddo !i enddo !j !$OMP END PARALLEL DO ! if (lpipe .and. lpipe_momtum .and. mslprf) then ! --- 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_sym1(util4, ip,text) util4(1:ii,1:jj) = mslprs(1:ii,1:jj,l1) write (text,'(a9,i3)') 'mslprs l=',l1 call pipe_compare_sym1(util4, ip,text) ! 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_sym1(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_sym1(util4, iv,text) ! 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_sym1(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_sym1(util4, iu,text) endif ! ! call dpudpv(dpu(1-nbdy,1-nbdy,1,m), ! & dpv(1-nbdy,1-nbdy,1,m), ! & p,depthu,depthv, max(0,margin-1)) !CL/RB MPI bug correction 2011-01 ! 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) ! ! --- account for temporal smoothing of mid-time dpmixl. calculate the vertical ! --- excursions of the coordinates immediately above and below the mixed ! --- layer base, then vertically interpolate this motion to dpmixl(i,j,m) ! if(hybrid .and. mxlkta) then ! ! --- rhs: dp, dpmixl.m ! --- lhs: util1, util2 ! margin = mbdy ! !$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 ! ! --- -------------- ! --- surface stress ! --- -------------- ! if (windf) then margin =0 #if defined (USE_NUOPC_CESMBETA) !jc weights for the coupled forcing if(cpl_implicit) then if(nstep2_cpl.eq.0 .or. (nstep1_cpl.eq.nstep2_cpl)) then cpl_w2=1. cpl_w3=0. else cpl_w2=(nstep-nstep1_cpl)/(nstep2_cpl-nstep1_cpl) cpl_w3=1.-cpl_w2 endif else cpl_w2=1.0 cpl_w3=0. endif #endif ! !$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 ! --- 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 ) ) 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)) 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 ! if (wndflg.eq.2 .or. wndflg.eq.3) then ! tau on p grid #if defined (USE_NUOPC_CESMBETA) if (cpl_taux) then surtx(i,j) = imp_taux(i,j,1)*cpl_w2 & + imp_taux(i,j,2)*cpl_w3 elseif (natm.eq.2) then surtx(i,j) = taux(i,j,l0)*w0+taux(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 endif ! cpl_taux if (cpl_tauy) then surty(i,j) = imp_tauy(i,j,1)*cpl_w2 & + imp_tauy(i,j,2)*cpl_w3 elseif (natm.eq.2) then surty(i,j) = tauy(i,j,l0)*w0+tauy(i,j,l1)*w1 else surty(i,j) = tauy(i,j,l0)*w0+tauy(i,j,l1)*w1 & + tauy(i,j,l2)*w2+tauy(i,j,l3)*w3 endif ! cpl_tauy #else 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 #endif /* USE_NUOPC_CESMBETA */ elseif (wndflg.eq.1) then ! tau on u&v grids - NOT RECOMMEDED #if defined (USE_NUOPC_CESMBETA) if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in momtum - wndflg should be ' write(lp,*) ' eq. to 2 with NUOPC' endif !1st tile call xcstop('(momtum)') stop '(momtum)' #else 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 #endif /* USE_NUOPC_CESMBETA */ else !wndflg.eq.4,5,6 ! --- calculate stress from 10m winds using cd_coare or cd_core2 ! --- for cd_core2, vpmx (vapmix) is specific humidity #if defined (USE_NUOPC_CESMBETA) if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in momtum - wndflg should be ' write(lp,*) ' eq. to 2 with NUOPC' endif !1st tile call xcstop('(momtum)') stop '(momtum)' #else 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 ! if (wndflg.eq.4 .and. flxflg.eq.6) then ! --- 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)) ! --- 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 ! --- vpmx assumed to contain specific humidity cdw = 1.0e-3*cd_core2(wind,vpmx,airt, & temp(i,j,1,n)) ! --- use U10 magnitude and direction for stress surtx(i,j) = rair*cdw*wind*wndx surty(i,j) = rair*cdw*wind*wndy !*****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 /* USE_NUOPC_CESMBETA */ endif !wndflg ! if (stroff) then surtx(i,j) = surtx(i,j) + oftaux(i,j) surty(i,j) = surty(i,j) + oftauy(i,j) endif ! 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 ! --- allow for ice-ocean stress #if defined (USE_NUOPC_CESMBETA) ! --- ice-ocean stress already in surtx and surty #else 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 /* USE_NUOPC_CESMBETA */ endif !ice-ocean stress endif !ip enddo !i enddo !j !$OMP END PARALLEL DO ! call xctilr(surtx,1,1, 6,6, halo_pv) call xctilr(surty,1,1, 6,6, halo_pv) endif !windf ! return ! contains include 'internal_kappaf.h' end subroutine momtum_hs ! real function cd_coare(wind,vpmx,airt,sst) implicit none ! real wind,vpmx,airt,sst ! ! --- Wind stress drag coefficient * 10^3 from an approximation ! --- to the COARE 3.0 bulk algorithm (Fairall et al. 2003). ! ! --- wind = wind speed (m/s) ! --- vpmx = water vapor mixing ratio (kg/kg) ! --- airt = air temperature (C) ! --- sst = sea temperature (C) ! ! --- Ta-Ts ! --- ============== ! --- STABLE: 7 to 0.75 degC ! --- NEUTRAL: 0.75 to -0.75 degC ! --- UNSTABLE: -0.75 to -8 degC ! ! --- Va ! --- ============== ! --- Low: 1 to 5 m/s ! --- High: 5 to 34 m/s ! ! --- vamax of 34 m/s from Sraj et al, 2013 (MWR-D-12-00228.1). ! real tamts,q,qva,va ! real, parameter :: vamin= 1.0, vamax=34.0 real, parameter :: tdmin= -8.0, tdmax= 7.0 real, parameter :: tzero=273.16 ! 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 ! ! --- 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)))))) ! ! --- 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 ! end function cd_coare ! real function cd_coarep(samo,vpmx,airt,pair,sst) implicit none ! real samo,vpmx,airt,pair,sst ! ! --- Wind stress drag coefficient * 10^3 from an approximation ! --- to the COARE 3.0 bulk algorithm (Fairall et al. 2003). ! ! --- samo = wind-ocean speed (m/s) ! --- vpmx = water vapor mixing ratio (kg/kg) ! --- airt = air temperature (C) ! --- pair = air pressure (Pa) ! --- sst = sea temperature (C) ! ! --- Ta-Ts ! --- ============== ! --- STABLE: 7 to 0.75 degC ! --- NEUTRAL: 0.75 to -0.75 degC ! --- UNSTABLE: -0.75 to -8 degC ! ! --- Va ! --- ============== ! --- Low: 1 to 5 m/s ! --- High: 5 to 34 m/s ! ! --- vamax of 34 m/s from Sraj et al, 2013 (MWR-D-12-00228.1). ! real tamts,q,qva,va ! real, parameter :: vamin= 1.0, vamax=34.0 real, parameter :: tdmin= -8.0, tdmax= 7.0 real, parameter :: tzero=273.16 ! 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 ! real satvpr,qsaturp,t,t6,p6 ! ! --- saturation vapor pressure (Pa), ! --- 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)))))) ! ! --- pressure dependent saturation mixing ratio (kg/kg) ! --- p6 is pressure in Pa qsaturp(t6,p6)=0.622*(satvpr(t6)/(p6-satvpr(t6))) ! ! --- 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 ! end function cd_coarep ! real function cd_core2(wind,sphm,airt,sst) implicit none ! real wind,sphm,airt,sst ! ! --- Wind stress drag coefficient * 10^3 from ! --- the CORE v2 bulk algorithm (Large and Yeager, 2009). ! ! --- wind = wind speed (m/s) ! --- sphm = specific humidity (kg/kg) ! --- airt = air temperature (C) ! --- sst = sea temperature (C) ! ! --- Added to HYCOM by Alexandra Bozec, FSU. ! 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 ! 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 ! ! --- saturation specific humidity real qsatur5,t,qra qsatur5(t,qra)= 0.98*qra*6.40380e5*exp(-5107.4/(t+tzero)) ! ! --- CORE v2 Large and Yeager 2009 Clim. Dyn.: The global climatology ! --- of an interannually varying air-sea flux dataset. ! --- The bulk formulae effectively transform the problem of specifying ! --- the turbulent surface fluxes (at zi=10m) into one of describing ! --- the near surface atmospheric state (wind, temperature and humidity). ! ! write(6,'(a,1p4g14.5)') ! & ' wind,sphm,airt,sst =',wind,sphm,airt,sst ! rair = 1.22 qrair = 1.0/rair zi = 10.0 tv = (airt+tzero)*(1.0+0.608*sphm) !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)* & (sphm-qsatur5(sst,qrair)) !L-Y eqn. 7c bstar = g*(tstar/tv+qstar/(sphm+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 ! 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 ! cd_core2 = cd10*1.0e3 return end function cd_core2 ! 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 implicit none ! integer m,n ! ! --- --------------------------------------------------------- ! --- momentum equations (2nd order version) ! ! --- Enstrophy conserving advection scheme (Sadourney, 1975) ! ! --- diffusion is Laplacian and/or biharmonic, both with ! --- "constant" and deformation dependent coefficients. ! ! --- hydrostatic equation and surface stress via momtum_hs ! --- --------------------------------------------------------- ! --- on entry: ! --- saln(:,:,:,m) = time step t with RA time smoothing ! --- saln(:,:,:,n) = time step t+1 ! ! --- dp(:,:,:,m) = time step t with RA time smoothing ! --- dp(:,:,:,n) = time step t+1 ! ! --- dpo(:,:,:,n) = time step t-1 ! --- dpo(:,:,:,m) = time step t ! ! --- dpv(:,:,:,n) = time step t-1 ! --- dpv(:,:,:,m) = time step t ! --- dpu(:,:,:,n) = time step t-1 ! --- dpu(:,:,:,m) = time step t ! ! --- v(:,:,:,n) = time step t-1 ! --- v(:,:,:,m) = time step t ! --- u(:,:,:,n) = time step t-1 ! --- u(:,:,:,m) = time step t ! ! --- on exit: ! --- dpv(:,:,:,m) = time step t with RA time smoothing ! --- dpv(:,:,:,n) = time step t+1 ! --- dpu(:,:,:,m) = time step t with RA time smoothing ! --- dpu(:,:,:,n) = time step t+1 ! ! --- v(:,:,:,m) = time step t with RA time smoothing ! --- v(:,:,:,n) = time step t+1 ! --- u(:,:,:,m) = time step t with RA time smoothing ! --- u(:,:,:,n) = time step t+1 ! ! --- vtotn(:,:) = time step t-1 to t+1 barotropic tendency ! --- utotn(:,:) = time step t-1 to t+1 barotropic tendency ! --- --------------------------------------------------------- ! logical, parameter :: lpipe_momtum=.false. !usually .false. ! #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 ! #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 ! #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,dzdx,dzdy ! real uatvk0,uatvkm,uatvkp,usd_v,vatuk0,vatukm,vatukp,vsd_u, & zbot,ztop #endif /* STOKES */ integer, save :: ifirst ! 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 real dp12,dp23,dp123,dp3m1,ql1,ql2,ql3 integer i,ia,ib,j,ja,jb,k,ka,l,mbdy,ktop,kmid,kbot,margin ! ! real*8 wtime ! external wtime ! real*8 wtime1(10),wtime2(20,kdm),wtimes ! character text*12 integer, save, allocatable, dimension(:,:) :: & mask ! real hfharm,a,b # include "stmt_fns.h" ! ! --- harmonic mean divided by 2 hfharm(a,b)=a*b/(a+b) ! data ifirst / 0 / ! #if defined(RELO) 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) ) call mem_stat_add( 2*(idm+2*nbdy)*(jdm+2*nbdy) ) usdflux = 0.0 !r_init vsdflux = 0.0 !r_init ! 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), & dzdx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk), & dzdy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,1:kk) ) call mem_stat_add( 6*(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 dzdx = 0.0 !r_init dzdy = 0.0 !r_init endif !usdflux #endif /* STOKES */ #endif /* RELO */ ! mbdy = 6 ! ! --- 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 (ifirst.eq.0) then ifirst=1 ! --- setup zero fill. margin = mbdy ! 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 ! dpthin = 0.001*onemm h1 = tenm !used in lateral weighting of hor.pres.grad. cutoff = 0.5*onem qcutoff = 1.0/cutoff ! ! --- --------------------------------------- ! --- hydrostatic equation and surface stress ! --- --------------------------------------- ! call momtum_hs(m,n) ! ! +++ ++++++++++++++++++ ! +++ momentum equations ! +++ ++++++++++++++++++ ! #if defined(STOKES) ! wtime1( 4) = wtime() margin = mbdy ! !$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 ! ! --- Calculate dudzu,dvdzu using a Lagrange interpolation. ! margin = mbdy-1 ! !$OMP PARALLEL DO PRIVATE(j,i,k, & !$OMP dp12,dp23,dp123,dp3m1,ql1,ql2,ql3, & !$OMP vatuk0,vatukm,vatukp, & !$OMP uatvk0,uatvkm,uatvkp) & !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin k=1 do i=1-margin,ii+margin if (SEA_U) then dp12 = max(cutoff,dpu(i,j,1,m))+ & max(cutoff,dpu(i,j,2,m)) dp23 = max(cutoff,dpu(i,j,2,m))+ & max(cutoff,dpu(i,j,3,m)) dp123 = dp12+dp23 ql1 = -2.0*(dp12+dp123)/(dp12 *dp123) ql2 = 2.0* dp123 /(dp12 *dp23) ql3 = -2.0* dp12 /(dp123*dp23) vatukm = 0.25*(v(i, j, 1,m)+ & v(i-1,j, 1,m)+ & v(i, j+1,1,m)+ & v(i-1,j+1,1,m) ) vatuk0 = 0.25*(v(i, j, 2, m)+ & v(i-1,j, 2, m)+ & v(i, j+1,2, m)+ & v(i-1,j+1,2, m) ) vatukp = 0.25*(v(i, j, 3,m)+ & v(i-1,j, 3,m)+ & v(i, j+1,3,m)+ & v(i-1,j+1,3,m) ) dudzu(i,j,1) = onem*(u(i,j,1,m)*ql1+ & u(i,j,2,m)*ql2+ & u(i,j,3,m)*ql3 ) dvdzu(i,j,k) = onem*(vatukm *ql1+ & vatuk0 *ql2+ & vatukp *ql3 ) endif !iu if (SEA_V) then dp12 = max(cutoff,dpv(i,j,1,m))+ & max(cutoff,dpv(i,j,2,m)) dp23 = max(cutoff,dpv(i,j,2,m))+ & max(cutoff,dpv(i,j,3,m)) dp123 = dp12+dp23 ql1 = -2.0*(dp12+dp123)/(dp12 *dp123) ql2 = 2.0* dp123 /(dp12 *dp23) ql3 = -2.0* dp12 /(dp123*dp23) uatvkm = 0.25*(u(i, j, 1,m)+ & u(i-1,j, 1,m)+ & u(i, j+1,1,m)+ & u(i-1,j+1,1,m) ) uatvk0 = 0.25*(u(i, j, 2, m)+ & u(i-1,j, 2, m)+ & u(i, j+1,2, m)+ & u(i-1,j+1,2, m) ) uatvkp = 0.25*(u(i, j, 3,m)+ & u(i-1,j, 3,m)+ & u(i, j+1,3,m)+ & u(i-1,j+1,3,m) ) dvdzv(i,j,1) = onem*(v(i,j,1,m)*ql1+ & v(i,j,2,m)*ql2+ & v(i,j,3,m)*ql3 ) dudzv(i,j,k) = onem*(uatvkm *ql1+ & uatvk0 *ql2+ & uatvkp *ql3 ) endif !iv enddo !i ! k=1 do k= 2,kk-1 do i=1-margin,ii+margin if (SEA_U) then dp12 = max(cutoff,dpu(i,j,k-1,m))+ & max(cutoff,dpu(i,j,k, m)) dp23 = max(cutoff,dpu(i,j,k, m))+ & max(cutoff,dpu(i,j,k+1,m)) dp123 = dp12+dp23 dp3m1 = dp23-dp12 ql1 = -2.0*dp23/(dp12*dp123) ql2 = 2.0*dp3m1/(dp12*dp23) ql3 = 2.0*dp12/(dp123*dp23) vatukm = 0.25*(v(i, j, k-1,m)+ & v(i-1,j, k-1,m)+ & v(i, j+1,k-1,m)+ & v(i-1,j+1,k-1,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, k+1,m)+ & v(i-1,j, k+1,m)+ & v(i, j+1,k+1,m)+ & v(i-1,j+1,k+1,m) ) dudzu(i,j,k) = onem*(u(i,j,k-1,m)*ql1+ & u(i,j,k, m)*ql2+ & u(i,j,k+1,m)*ql3 ) dvdzu(i,j,k) = onem*(vatukm* ql1+ & vatuk0 *ql2+ & vatukp *ql3 ) !diag if (max(abs(i-itest),abs(j-jtest)).eq.0) then !diag write(lp,'(a,i3,3g20.10)') & !diag 'dp,k =',k,dp12,dp23,ql1+ql2+ql3 !diag write(lp,'(a,i3,3g20.10)') & !diag 'ql,k =',k,ql1,ql2,ql3 !diag write(lp,'(a,i3,2g20.10)') & !diag 'dz,k =',k,dudzu(i,j,k),dvdzu(i,j,k) !diag endif !test endif !iu if (SEA_V) then dp12 = max(cutoff,dpv(i,j,k-1,m))+ & max(cutoff,dpv(i,j,k, m)) dp23 = max(cutoff,dpv(i,j,k, m))+ & max(cutoff,dpv(i,j,k+1,m)) dp123 = dp12+dp23 dp3m1 = dp23-dp12 ql1 = -2.0*dp23/(dp12*dp123) ql2 = 2.0*dp3m1/(dp12*dp23) ql3 = 2.0*dp12/(dp123*dp23) uatvkm = 0.25*(u(i, j, k-1,m)+ & u(i-1,j, k-1,m)+ & u(i, j+1,k-1,m)+ & u(i-1,j+1,k-1,m) ) uatvk0 = 0.25*(u(i, j, k, m)+ & u(i-1,j, k, m)+ & u(i, j+1,k, m)+ & u(i-1,j+1,k, m) ) uatvkp = 0.25*(u(i, j, k+1,m)+ & u(i-1,j, k+1,m)+ & u(i, j+1,k+1,m)+ & u(i-1,j+1,k+1,m) ) dvdzv(i,j,k) = onem*(v(i,j,k-1,m)*ql1+ & v(i,j,k ,m)*ql2+ & v(i,j,k+1,m)*ql3 ) dudzv(i,j,k) = onem*(uatvkm *ql1+ & uatvk0 *ql2+ & uatvkp *ql3 ) endif !iv enddo !i enddo !k k=kk do i=1-margin,ii+margin if (SEA_U) then dp12 = max(cutoff,dpu(i,j,kk-2,m))+ & max(cutoff,dpu(i,j,kk-1,m)) dp23 = max(cutoff,dpu(i,j,kk-1,m))+ & max(cutoff,dpu(i,j,kk, m)) dp123 = dp12+dp23 ql1 = 2.0*dp23/(dp12*dp123) ql2 = -2.0*dp123/(dp12*dp23) ql3 = 2.0*(dp12+dp123)/(dp123*dp23) vatukm = 0.25*(v(i, j ,kk-2,m)+ & v(i-1,j ,kk-2,m)+ & v(i, j+1,kk-2,m)+ & v(i-1,j+1,kk-2,m) ) vatuk0 = 0.25*(v(i, j ,kk-1, m)+ & v(i-1,j ,kk-1, m)+ & v(i, j+1,kk-1, m)+ & v(i-1,j+1,kk-1, m) ) vatukp = 0.25*(v(i, j, kk,m)+ & v(i-1,j, kk,m)+ & v(i, j+1,kk,m)+ & v(i-1,j+1,kk,m) ) dudzu(i,j,kk) = onem*(u(i,j,kk-2,m)*ql1+ & u(i,j,kk-1,m)*ql2+ & u(i,j,kk, m)*ql3 ) dvdzu(i,j,kk) = onem*(vatukm *ql1+ & vatuk0 *ql2+ & vatukp *ql3 ) endif !iu if (SEA_V) then dp12 = max(cutoff,dpv(i,j,kk-2,m))+ & max(cutoff,dpv(i,j,kk-1,m)) dp23 = max(cutoff,dpv(i,j,kk-1,m))+ & max(cutoff,dpv(i,j,kk, m)) dp123 = dp12+dp23 ql1 = 2.0*dp23/(dp12*dp123) ql2 = -2.0*dp123/(dp12*dp23) ql3 = 2.0*(dp12+dp123)/(dp123*dp23) uatvkm = 0.25*(u(i, j ,kk-2,m)+ & u(i-1,j ,kk-2,m)+ & u(i, j+1,kk-2,m)+ & u(i-1,j+1,kk-2,m) ) uatvk0 = 0.25*(u(i, j ,kk-1,m)+ & u(i-1,j ,kk-1,m)+ & u(i, j+1,kk-1,m)+ & u(i-1,j+1,kk-1,m) ) uatvkp = 0.25*(u(i, j, kk, m)+ & u(i-1,j, kk, m)+ & u(i, j+1,kk, m)+ & u(i-1,j+1,kk, m) ) dvdzv(i,j,kk) = onem*(v(i,j,kk-2,m)*ql1+ & v(i,j,kk-1,m)*ql2+ & v(i,j,kk, m)*ql3 ) dudzv(i,j,kk) = onem*(uatvkm *ql1+ & uatvk0 *ql2+ & uatvkp *ql3 ) endif !iv enddo !i ! k=kk enddo !j ! ! Calculation of dzdx,dzdy needed in momtum -u-v- ! ----------------------------------------- ! margin = 0 ! !$OMP PARALLEL DO PRIVATE(j,i,k,zbot,ztop) & !$OMP SCHEDULE(STATIC,jblk) do j = 1-margin, jj+margin do i=1-margin,ii+margin if (SEA_U) then ztop = -(depths(i,j)-depths(i-1,j)) do k = kk, 1, -1 zbot = ztop ztop = zbot + & qonem*(dp(i, j,k,m)*oneta(i, j,m)- & dp(i-1,j,k,m)*oneta(i-1,j,m) ) dzdx(i,j,k) = 0.5*(ztop+zbot) enddo !k endif ! iu ! if (SEA_V) then ztop = -(depths(i,j)-depths(i,j-1)) do k = kk, 1, -1 zbot = ztop ztop = zbot + & qonem*(dp(i,j, k,m)*oneta(i, j,m)- & dp(i,j-1,k,m)*oneta(i,j-1,m) ) dzdy(i,j,k) = 0.5*(ztop+zbot) enddo !k endif ! iv enddo !i enddo !j 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 */ ! ! --- rhs: p, u.n+, v.n+, ubavg.n+, vbavg.n+, depthv+, pvtrop+ ! --- rhs: dpmixl.m+, taux+, dpu, depthu+, dpv, tauy+ ! --- lhs: util1, util2, drag, ubrhs, stresx, vbrhs, stresy ! !-----disp_count = disp_count + 1 stresl(:,:) = 0.0 ! if (drglim.gt.0.0) then adrlim = drglim else adrlim = 0.125 endif dt1inv = 1./delt1 ! margin = mbdy - 1 ! !$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 ! do i=1-margin,ii+margin if (SEA_P) then ! ! --- bottom drag (standard bulk formula) ! --- bottom stress is applied over thickness thkbot (or the ! --- total depth is this is less) ! thkbop(i,j)=min( thkbot*onem, p(i,j,kk+1) ) ! ! --- the bottom stress term is estimated using velocity averaged ! --- over the bottom boundary layer. this thickness is dpbbl for ! --- the kpp boundary layer; otherwise, it is thkbop ! 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)) ! ! --- drag = Cb * (|v| + c.bar) ! --- include 1/thkbop for the fraction of layer calculation ! --- and onem for conversion from 1/dp to 1/h ! vmag=0.5*sqrt(ubot**2+vbot**2) !0.5 from u+u.i-1 & v+v.j-1 dall=cbp(i,j)*(vmag+cbarp(i,j)) !no tidal drag drag(i,j)=dall*onem/thkbop(i,j) util1(i,j) = ubot util2(i,j) = vbot if (mxlkpp .and. bblkpp) then ustarb(i,j)=sqrt(dall*vmag) endif ! ! --- tidal bottom drag, drgten.1.1 is drgscl*rh ! 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 ! ! --- store r.h.s. of barotropic u/v eqn. in -ubrhs,vbrhs- ! --- time-interpolate wind stress ! 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 ! 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 ! --- units of surtx are N/m^2 (i.e. Pa) stresx(i,j)=(surtx(i,j)+surtx(i-1,j))*0.5*g & /(pstres*svref) else ! no taux stresx(i,j)=0. endif !windf:else endif !iu ! 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 ! 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 ! --- units of surty are N/m^2 (i.e. Pa) stresy(i,j)=(surty(i,j)+surty(i,j-1))*0.5*g & /(pstres*svref) else ! no tauy stresy(i,j)=0. endif !windf:else endif !iv enddo !i enddo !j !$OMP END PARALLEL DO ! if (lpipe .and. lpipe_momtum) then ! --- compare two model runs. write (text,'(a9,i3)') 'uba.n n=',n call pipe_compare_sym1(ubavg(1-nbdy,1-nbdy,n),iu,text) write (text,'(a9,i3)') 'vba.n n=',n call pipe_compare_sym1(vbavg(1-nbdy,1-nbdy,n),iv,text) write (text,'(a9,i3)') 'cbarp k=',0 call pipe_compare_sym1(cbarp,ip,text) write (text,'(a9,i3)') 'cbp k=',0 call pipe_compare_sym1(cbp ,ip,text) write (text,'(a9,i3)') 'thkbop k=',0 call pipe_compare_sym1(thkbop,ip,text) write (text,'(a9,i3)') 'ubot_p k=',0 call pipe_compare_sym1(util1 ,ip,text) write (text,'(a9,i3)') 'vbot_p k=',0 call pipe_compare_sym1(util2 ,ip,text) write (text,'(a9,i3)') 'drag k=',0 call pipe_compare_sym1(drag,ip,text) endif ! ! --- the old momeq2.f starts here ! ! wtime1( 5) = wtime() ! ! --- rhs: 0.0 ! --- lhs: util1, util2 ! margin = mbdy - 2 ! !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin ! --- spatial weighting function for pressure gradient calculation: util1(i,j)=0.0 util2(i,j)=0.0 ! --- p.1 pnkp( i,j)=0.0 enddo !i enddo !j ! dpmx(:,:) = 2.*cutoff !otherwise ii(jj)+1 to i(j)dm are NaN ! do 9 k=1,kk ! ! --- store total (barotropic plus baroclinic) flow at old and mid time in ! --- -utotn,vtotn- and -utotm,vtotm- respectively. store minimum thickness ! --- values for use in pot.vort. calculation in -dpmx-. ! --- store p.k and p.k+1 for time level t+1 in pnk0,pnkp ! ! wtime2( 1,k) = wtime() ! ! --- rhs: dpmx, dp.m+ ! --- lhs: dpmx ! margin = mbdy - 2 ! 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 ! ! wtime2( 2,k) = wtime() ! ! --- rhs: ubavg.m, ubavg.n, dp.m+, dpu ! --- lhs: utotm, utotn, uflux, dpmx, pu ! margin = mbdy - 2 ! !$OMP PARALLEL DO PRIVATE(j,l,i) & !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do l=1,isu(j) !ok ! 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)* & 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)* & max(dpo(i-1,j,k,m),cutoff) #endif endif endif ! 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)*max(dpu(i,j,k,m),cutoff) !usd is total drift #endif pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,m) ! enddo !i enddo !l enddo !j !$OMP END PARALLEL DO ! ! wtime2( 3,k) = wtime() ! ! --- rhs: vbavg.m, vbavg.n, dp.m+, dpv ! --- lhs: vtotm, vtotn, vflux, dpmx, pv ! margin = mbdy - 2 ! 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)* & 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)* & max(dpo(i,j-1,k,m),cutoff) #endif endif endif enddo !l enddo !i ! ! 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)*max(dpv(i,j,k,m),cutoff) !vsd is total drift #endif pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,m) endif !iv enddo !i enddo !j ! ! --- define auxiliary velocity fields (via,vib,uja,ujb) to implement ! --- sidewall friction along near-vertical bottom slopes. wgtja,wgtjb,wgtia, ! --- wgtib indicate the extent to which a sidewall is present. ! ! wtime2( 5,k) = wtime() ! ! --- rhs: pu, depthu+, utotn+, wgtja ! --- lhs: wgtja, wgtjb, uja, ujb, dl2u ! --- rhs: pv, depthv+, vtotn+, wgtia ! --- lhs: wgtia, wgtib, via, vib, dl2v ! --- rhs: vtotm, vort+, corio+, dp.m+, dpmx+, vtotn ! --- lhs: vort, potvor, defor2 ! margin = mbdy - 2 ! !$OMP PARALLEL DO PRIVATE(j,ja,jb,l,i,ia,ib,aspy2,aspx2) & !$OMP SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin ! --- assume margin to properly commute pressure torque from layer to layer !> Jul. 2000 - loop reordering and logic changes for OpenMP !> Aug. 2000 - loop 117 executed only when hybrid vertical coordinate is used !> Oct. 2001 - replaced biharm with veldf[24] and visco[24] !> Sep. 2004 - kapref selects one of three thermobaric reference states !> Jun. 2006 - split out momtum_hs, for initial ssh calculation !> Nov. 2006 - inserted volume-force for tide !> Apr. 2007 - added drglim, implicit or CFL-limited explicit bottom drag !> Apr. 2007 - added dragrh, linear tidal drag based on bottom roughness !> Apr. 2007 - btrlfr: moved [uvp]bavg assignment to barotp !> Jun. 2007 - added momtum4 !> Mar 2009 - more accurate kappaf, with potential density !> Now 2009 - several bugfixes to momtum4 !> Mar 2010 - remove anti-drag on detided velocity !> Apr 2010 - put back anti-drag on detided velocity !> Apr 2011 - cbarp(i,j) in place of cbar !> Jul 2011 - salfac(i,j) in place of tidsal !> Jul 2011 - modified momtum4 based on NERSC 2.1.03_MPI momtum_quick !> Jul 2011 - modified momtum4 to use viscp and viscq !> Jul 2011 - modified momtum4 to use local arrays for y advection !> Jul 2011 - modified momtum4 to always use v 0 -v for extrapolation !> Aug 2011 - use ra2fac for wuv1 and wuv2 (so now wuv==wts) !> Aug 2011 - replaced dpold,dpoldm with dpo !> Aug 2011 - reworked Robert-Asselin time filter !> Sep 2011 - cbp(i,j) in place of cb !> Jan 2012 - added thkcdw !> July 2012 - bugfix for bottom drag when depth < thkbot !> July 2012 - thkbop is now always based on thkbot (even for bblkpp) !> Nov. 2012 - surtx,y are halo_pv (not halo_ps) !> Nov. 2012 - wndflg=4 for calculating wind stress using cd_coare !> Nov. 2012 - oftaux,oftauy for wind stress offset !> Jan. 2013 - added thkdrg=0.0 for tidal drag in barotp.f !> Jan. 2013 - replaced dragrh with drgten.1.1 !> Jan. 2013 - added gtide for tidal body forcing !> May 2013 - removed thbase from montg.kk !> July 2013 - vamax set to 34 m/s !> Sep. 2013 - optionally added Stokes drift !> Nov. 2013 - wndflg=5 for calculating wind stress using cd_core2 !> Jan. 2014 - tv is in Kelvin (cd_core2) !> Jan. 2014 - added gslpr for atmospheric pressure forcing !> Jan. 2014 - added natm !> Apr. 2014 - added pair for atmospheric pressure forcing !> Apr. 2014 - added ice shelf logic (ishlf) !> May 2014 - use land/sea masks (e.g. ip) to skip land !> Oct 2014 - added cd_coarep for flxflg=6 !> Oct 2014 - flxflg=6 uses samo in place of wind everywhere !> Oct 2014 - flxflg=4,5 uses wind-current magnitude and direction for stress !> Oct. 2042 - dragw_rho appropriate for thkcdw=3.0, 2x default !> Apr. 2015 - moved atmospheric pressure forcing to barotp !> Apr. 2015 - moved tidal body forcing to barotp !> Sep. 2015 - calculate dudzu,dvdzu using a Lagrange interpolation !> July 2017 - backported momtum4 from 2.2.101 - start comments !> Mar. 2016 - momtum4 diffuson uses lagged layer thickness !> Mar. 2016 - momtum4 uses coefficient arrays to handle near-land stencils !> Mar. 2016 - addded momtyp==3 to momtum4 for Split QUICK advection !> Mar. 2016 - montum4 handles non-constant grids !> Mar. 2016 - montum4 allows non-zero visco2 and veldf2 !> Mar. 2016 - montum4 allows non-zero visco2 and veldf2 !> Mar. 2016 - in montum4, veldf[24]u is actually veldf[24]p !> July 2017 - backported momtum4 from 2.2.101 - end comments !> July 2017 - added momtum4_cfl compile time flag !> July 2017 - bugfix for qrair in cd_core2 !> Aug 2017 - kappaf via internal function !> Apr 2018 - cd_core2 uses specific humidity, not mixing ratio !> Apr 2018 - wndflg=5 uses U10, not U10-Uocn !> Oct 2018 - centered pressure for kapref available via a compile-time macro !> Dec 2018 - add /* USE_NUOPC_CESMBETA */ macro for coupled simulation !> Dec 2018 - use max (not average) to map drag from p to u and v grids !> Feb. 2019 - added montg_c correction to psikk (sshflg.eq.2) !> Sep. 2019 - added momtum_init