#ifdef HAVE_CONFIG_H #include "config.h" #endif module shallow_water_mod ! ------------------------ use kinds, only : real_kind ! ------------------------ use physical_constants, only : omega, g, rearth=>rearth0, rrearth=>rrearth0, dd_pi use physical_constants, only : Lx, Ly ! ------------------------ use dimensions_mod, only : nlev, np ! ------------------------ use derivative_mod, only : derivative_t, vorticity_sphere, & divergence, vorticity, laplace_sphere_wk, & vlaplace_sphere_wk, divergence_sphere ! ------------------------ use edge_mod, only : edgeVpack, edgeVunpack use edgetype_mod, only : EdgeBuffer_t ! ------------------------ use bndry_mod, only : bndry_exchangeV ! ------------------------ use coordinate_systems_mod, only : spherical_polar_t ! ------------------------ use quadrature_mod, only : gauss, gaussian_int, quadrature_t ! ------------------------ use global_norms_mod, only : global_integral, linf_snorm, l1_snorm, l2_snorm ! ------------------------ use time_mod, only : time_at, secpday, timelevel_t, nmax ! ------------------------ use element_mod, only : element_t ! ------------------------ use hybrid_mod, only : hybrid_t ! ------------------------ use parallel_mod, only : parallel_t ! ------------------------ use common_io_mod, only: output_prefix ! Added to support output_prefix ! ------------------------ use viscosity_mod, only: biharmonic_wk, neighbor_minmax ! ------------------------ use control_mod, only : nu, nu_div, nu_s, hypervis_order, hypervis_subcycle, limiter_option, integration, test_case, sub_case, & kmass, g_sw_output,toy_chemistry ! ------------------------ implicit none private logical :: compute_penst = .true. ! ============================================ ! Shallow Water Test Case Parameters ! for test case 1 and 2; Cosine Bell and ! Steady State Geostrophic ! ============================================ real (kind=real_kind), private, parameter :: h0 = 1000.0D0 ! height of cosine bell in meters real (kind=real_kind), private, parameter :: gh0_tc1= g*h0 ! geopotential height of cosine bell real (kind=real_kind), private, parameter :: rr = 1.0/3.0 ! radius of cosine bell real (kind=real_kind), public, parameter :: latc0 = 0.0D0 ! init latitude of cosine bell center real (kind=real_kind), public, parameter :: lonc0 = 1.5D0*DD_PI ! init longitude of cosine bell center real (kind=real_kind), private, parameter :: u0 = 2.0D0*DD_PI*rearth/(12.0*86400.0) ! adv. wind velocity real (kind=real_kind), private, parameter :: alpha = 0.0D0 ! angle of wind advection rel. to equator ! real (kind=real_kind), private, parameter :: alpha = DD_PI/2.0D0 ! angle of wind advection rel. to equator ! real (kind=real_kind), private, parameter :: alpha = DD_PI/4.0D0 ! angle of wind advection rel. to equator real (kind=real_kind), private, parameter :: gh0_tc2 = 2.94D4 ! m^2/s^2 (p. 218 of Williamson) ! ================================================ ! Shallow Water Test Case Parameters ! Parameters for test case 5, Flow over Mountain ! ================================================ real (kind=real_kind), private, parameter :: h0_tc5 = 5960.0D0 ! height field real (kind=real_kind), parameter, private :: gh0_tc5 = g*h0_tc5 ! test case 5 geopotential real (kind=real_kind), private, parameter :: u0_tc5 = 20.0D0 ! wind velocity real (kind=real_kind), private, parameter :: h_mtn = 2000.0D0 ! height of mountain real (kind=real_kind), private, parameter :: R_mtn = DD_PI/9.0D0 ! radius of mountain real (kind=real_kind), public, parameter :: lat_mtn = DD_PI/6.0D0 ! latitude of mountain real (kind=real_kind), public, parameter :: lon_mtn = 1.5D0*DD_PI ! longitude of mountain ! ============================================ ! Shallow Water Test Case Parameters ! Parameters for test case 6, Rossby-Haurwitz ! ============================================ real (kind=real_kind), parameter, private :: h0_tc6 = 8.0D3 real (kind=real_kind), parameter, private :: gh0_tc6 = g*h0_tc6 real (kind=real_kind), parameter, private :: omega_rh = 7.848D-6 real (kind=real_kind), parameter, private :: K_rh = omega_rh real (kind=real_kind), private :: R_rh = 4 ! ============================================ ! Shallow Water Test Case Parameters ! Parameters for test case 8, ! Barotropic Instability. ! ============================================ real (kind=real_kind), parameter, private :: h0_tc8 = 10000.0D0 ! mean fluid depth (meters) real (kind=real_kind), parameter, private :: hhat_tc8 = 120.0D0 ! perturbation depth (meters) real (kind=real_kind), parameter, private :: gh0_tc8 = g*h0_tc8 ! mean geopotential real (kind=real_kind), private, parameter :: u0_tc8 = 80.0D0 ! wind velocity real (kind=real_kind), private, parameter :: alpha_tc8 = 3.0D0 ! meridional scale factor in perturbation real (kind=real_kind), private, parameter :: beta_tc8 = 15.0D0 ! zonal scale factor in perturbation real (kind=real_kind), public, parameter :: lat_tc8 = DD_PI/4.0D0 ! latitude of perturbation real (kind=real_kind), public, parameter :: lon_tc8 = DD_PI/4.0D0 ! longitude of perturbation real (kind=real_kind), public, parameter :: gamma_tc8 = DD_PI/18.0D0 ! zonal velocity field scale factor real (kind=real_kind), public, parameter :: eps_tc8 = 1.0D-10 ! relative error of balance integral integer, parameter :: Ngs_tc8 = 128 ! number of gauss points in balance integral ! ========================================== ! Vortex problem: Nair & Machenhaur, 2002 ! ========================================== real (kind=real_kind), private, parameter :: alpha_vortex =DD_PI*0.25D0 ! meridional scale factor in perturbation real (kind=real_kind), parameter, private :: a_vortex = rearth real (kind=real_kind), private, parameter :: u0_vortex = 2*DD_PI*a_vortex/10.368D5 ! wind velocity real (kind=real_kind), private, parameter :: omega_s_vortex = u0_vortex/a_vortex real (kind=real_kind), private, parameter :: lat0_vort1_vortex=0; real (kind=real_kind), private, parameter :: lon0_vort1_vortex=DD_PI*0.5D0; !Vortex center at time t=0 logical, private, parameter :: static_vortex = .false. ! ============================================ ! Shallow Water Test Case Parameters ! Parameters for test case SJ1, ! Strong jet northern hemisphere. ! ============================================ real (kind=real_kind), parameter, private :: h0_sj1 = 10000.0D0 ! mean fluid depth (meters) real (kind=real_kind), parameter, private :: hhat_sj1 = 120.0D0 ! perturbation depth (meters) real (kind=real_kind), private :: gh0_sj1 = g*h0_sj1 ! mean geopotential real (kind=real_kind), private, parameter :: u0_sj1 = 80.0D0 ! wind velocity real (kind=real_kind), private, parameter :: alpha_sj1 = 3.0D0 ! meridional scale factor in perturbation real (kind=real_kind), private, parameter :: beta_sj1 = 15.0D0 ! zonal scale factor in perturbation real (kind=real_kind), public, parameter :: lat_sj1 = DD_PI/4.0D0 ! latitude of perturbation real (kind=real_kind), public, parameter :: lon_sj1 = DD_PI/4.0D0 ! longitude of perturbation real (kind=real_kind), public, parameter :: gamma_sj1 = DD_PI/18.0D0 ! zonal velocity field scale factor real (kind=real_kind), public, parameter :: eps_sj1 = 1.0D-10 ! relative error of balance integral integer, parameter :: Ngs_sj1 = 200 ! number of gauss points in balance integral real (kind=real_kind), public, parameter :: lat0_sj1 = DD_PI/7.0D0 ! used for velocity real (kind=real_kind), public, parameter :: lat1_sj1 = 5.0D0*DD_PI/14.0D0! used for velocity real (kind=real_kind), public, parameter :: lat2_sj1 = DD_PI/4.0D0 ! used for velocity ! ========================================== ! Swirl problem: Nair & Lauritzen, 2010 ! ========================================== integer, public, parameter :: kmass_swirl=4 !at which level the const field sits real (kind=real_kind), private, parameter :: mean_pressure=0.0D0 real (kind=real_kind), private, parameter :: bellradius=0.50D0 !gauss bell would use different radius real (kind=real_kind), private, parameter :: rotangle=0.0D0 ! IN RADIANS real (kind=real_kind), private, parameter :: tracer_lowest=0.1D0 real (kind=real_kind), private, parameter :: tracer_highest=0.9D0 logical :: add_pure_rotation real (kind=real_kind) :: lat1, lon1, lat2, lon2, Tperiod, Kcoef !sub_case = 1 real (kind=real_kind), private, parameter :: lon1_case1=DD_PI real (kind=real_kind), private, parameter :: lat1_case1=DD_PI/3.0D0; real (kind=real_kind), private, parameter :: lon2_case1=DD_PI real (kind=real_kind), private, parameter :: lat2_case1=-DD_PI/3.0D0; !sub_case = 2 real (kind=real_kind), private, parameter :: lon1_case2=5.0D0*DD_PI/6.0D0; real (kind=real_kind), private, parameter :: lat1_case2=0.0D0; real (kind=real_kind), private, parameter :: lon2_case2=7.0D0*DD_PI/6.0D0; real (kind=real_kind), private, parameter :: lat2_case2=0.0D0; !sub_case = 3 !real (kind=real_kind), private, parameter :: lon1_case3=3.0D0*DD_PI/4.0D0 !real (kind=real_kind), private, parameter :: lat1_case3=0.0D0; !real (kind=real_kind), private, parameter :: lon2_case3=5.0D0*DD_PI/4.0D0 !real (kind=real_kind), private, parameter :: lat2_case3=0.0D0; ! i put subcase 2 center data here cause this is what pl uses for paper real (kind=real_kind), private, parameter :: lon1_case3=5.0D0*DD_PI/6.0D0; real (kind=real_kind), private, parameter :: lat1_case3=0.0D0; real (kind=real_kind), private, parameter :: lon2_case3=7.0D0*DD_PI/6.0D0; real (kind=real_kind), private, parameter :: lat2_case3=0.0D0; !subcase 4 is the same as subcase 2 but with added zonal flow, see below ! ============================================ ! Shallow Water Test Case Parameters ! Parameters for test case planar_dbl_vrtx, ! Planar Double Vortex ! ============================================ real (kind=real_kind), parameter, private :: h0_dbv = 750.0D0 ! mean fluid depth (meters) real (kind=real_kind), parameter, private :: dh_dbv = 75.0D0 ! perturbation depth (meters) real (kind=real_kind), private :: gh0_dbv = g*h0_dbv ! mean geopotential real (kind=real_kind), private, parameter :: f_dbv = 0.00006147D0 ! Coriolis parameter real (kind=real_kind), private, parameter :: sigmax_dbv = 3.0D0/40.0D0 ! vortex size in y-dir as a fraction of Lx real (kind=real_kind), private, parameter :: sigmay_dbv = 3.0D0/40.0D0 ! vortex size in x-dir as a fraction of Ly real (kind=real_kind), private, parameter :: ox_dbv = 0.1D0 ! offset in x-dir for vortices real (kind=real_kind), private, parameter :: oy_dbv = 0.1D0 ! offset in y-dir for vortices public :: sweq_invariants public :: tc1_init_state ! Initialize test case 1: Cosine Bell public :: tc1_init_pmean ! Initialize pmean for test case 1 public :: tc1_velocity public :: tc1_geopotential public :: tc1_phi ! test case 1 on unstaggered grid (nair) public :: grad_cosine_bell public :: tc1_errors public :: tc2_init_state ! Initialize test case 2: Global steady state nonlinear geostrophic flow public :: tc2_init_pmean public :: tc2_geopotential private :: tc2_coriolis_init public :: tc2_errors public :: tc2_phi public :: tc5_init_state ! Initialize test case 5: Flow Over a Mountain public :: tc5_init_pmean public :: tc5_velocity private :: tc5_coriolis_init public :: tc5_mountain public :: tc5_geopotential public :: tc5_errors public :: tc5_invariants public :: tc5_phi #if 0 public :: tc5_init_diag #endif public :: tc6_init_state ! Initialize test case 6: Rossby-Haurwitz Wave public :: tc6_init_pmean private :: tc6_geopotential private :: tc6_velocity public :: tc6_errors public :: tc8_init_state ! Initialize test case 8: Barotropic Instability public :: tc8_init_pmean private :: tc8_balance private :: tc8_perturbation private :: tc8_integrand private :: tc8_velocity public :: vortex_init_state public :: vortex_velocity public :: vortex_errors public :: vortex_exact public :: swirl_velocity public :: swirl_init_state public :: swirl_errors public :: toy_chemistry_forcing public :: sj1_init_state ! Initialize strong jet case 1: Galewski et al. public :: sj1_init_pmean private :: sj1_balance private :: sj1_perturbation private :: sj1_integrand private :: sj1_velocity private :: sj1_velocity_cubedsphere public :: sj1_errors !public :: sj1_invariants public :: planar_dbl_vrtx_init_state ! Initialize planar double vortex test case private :: planar_dbl_vrtx_geopotential private :: planar_dbl_vrtx_velocity public :: planar_dbl_vrtx_invariants real (kind=real_kind), public,save :: Imass0,Ienergy0,Ipenst0 ! common /comintegral/Imass0,Ienergy0,Ipenst0 real (kind=real_kind), public,save :: Ibalance_tot ! common /comtc8/Ibalance_tot contains ! ====================================================== ! shallow water equation invariants on a non-staggered grid. ! ====================================================== subroutine sweq_invariants(elem, iounit,tl,pmean,edge3,deriv,hybrid,nets,nete) type(element_t), intent(inout), target :: elem(:) integer , intent(in) :: iounit type (TimeLevel_t) , intent(in) :: tl real (kind=real_kind), intent(in) :: pmean type (EdgeBuffer_t) :: edge3 type (derivative_t) :: deriv type (hybrid_t) , intent(in) :: hybrid integer , intent(in) :: nets,nete ! Local variables integer :: ie,i,j,kptr real (kind=real_kind) :: vco(np,np,2) real (kind=real_kind) :: gv(np,np,2) real (kind=real_kind) :: E(np,np) real (kind=real_kind) :: mass(np,np,nlev,nets:nete) real (kind=real_kind) :: kenergy(np,np,nets:nete) real (kind=real_kind) :: penergy(np,np,nets:nete) real (kind=real_kind) :: penst(np,np,nets:nete) real (kind=real_kind) :: pv(np,np,nets:nete) real (kind=real_kind) :: div(np,np,nets:nete) real (kind=real_kind) :: vor(np,np,nets:nete) real (kind=real_kind) :: diss_p(np,np,nlev,nets:nete) real (kind=real_kind) :: diss_v(np,np,2,nlev,nets:nete) real (kind=real_kind) :: v1,v2,vlatlon(np,np,2) real (kind=real_kind) :: h,hstar,hs real (kind=real_kind) :: dm,dKE,dPE,dENS,dissE ! dont use Imass0,Ienergy0,Ipenst0 because they are being used by tc5_invariants real (kind=real_kind),save :: Imass_init(nlev),Ienergy_init,Ipenst_init,Ipv_init,Idiv_init real (kind=real_kind),save :: time_last,Imass_last(nlev),Ikenergy_last,Ipenergy_last,Ipenst_last,Ipv_last,Idiv_last real (kind=real_kind) :: time,Imass(nlev),Ikenergy,Ipenergy,Ipenst,Ienergy,Ipv,Idiv,rms_div,rms_vor,rms_f real (kind=real_kind), dimension(:,:,:), pointer :: v integer :: k,n0,k1 n0 = tl%n0 k=1 if (tl%nstep == 0) then if (hybrid%masterthread) then open(iounit+0,file=TRIM(output_prefix)//"sweq.mass",status="unknown",form="formatted") open(iounit+1,file=TRIM(output_prefix)//"sweq.ke",status="unknown",form="formatted") open(iounit+2,file=TRIM(output_prefix)//"sweq.pe",status="unknown",form="formatted") open(iounit+3,file=TRIM(output_prefix)//"sweq.penst",status="unknown",form="formatted") open(iounit+4,file=TRIM(output_prefix)//"sweq.pv",status="unknown",form="formatted") open(iounit+5,file=TRIM(output_prefix)//"sweq.div",status="unknown",form="formatted") open(iounit+6,file=TRIM(output_prefix)//"sweq.div2",status="unknown",form="formatted") open(iounit+7,file=TRIM(output_prefix)//"sweq.vor2",status="unknown",form="formatted") endif endif if (nu > 0 .or. nu_s > 0 ) then ! convert to lat-lon (needed to call biharmonic() ) do ie=nets,nete v => elem(ie)%state%v(:,:,:,k,n0) do j=1,np do i=1,np v1 = elem(ie)%state%v(i,j,1,k,n0) ! contra v2 = elem(ie)%state%v(i,j,2,k,n0) ! contra v(i,j,1)=elem(ie)%D(i,j,1,1)*v1 + elem(ie)%D(i,j,1,2)*v2 ! contra->latlon v(i,j,2)=elem(ie)%D(i,j,2,1)*v1 + elem(ie)%D(i,j,2,2)*v2 ! contra->latlon enddo enddo enddo if (hypervis_order==2) then call biharmonic_wk(elem,diss_p,diss_v,deriv,edge3,hybrid,n0,nets,nete) else do ie=nets,nete v => elem(ie)%state%v(:,:,:,k,n0) E(:,:)=elem(ie)%state%p(:,:,k,n0) + elem(ie)%state%ps(:,:) diss_p(:,:,k,ie)=laplace_sphere_wk(E,deriv,elem(ie),var_coef=.false.) diss_v(:,:,:,k,ie)=vlaplace_sphere_wk(v,deriv,elem(ie),var_coef=.false.) enddo endif ! convert lat-lon -> contra variant do ie=nets,nete v => elem(ie)%state%v(:,:,:,k,n0) do j=1,np do i=1,np v1=v(i,j,1) v2=v(i,j,2) v(i,j,1) = elem(ie)%Dinv(i,j,1,1)*v1 + elem(ie)%Dinv(i,j,1,2)*v2 v(i,j,2) = elem(ie)%Dinv(i,j,2,1)*v1 + elem(ie)%Dinv(i,j,2,2)*v2 enddo enddo enddo else diss_p=0; diss_v=0; endif do ie=nets,nete do j=1,np do i=1,np k1=1 v1 = elem(ie)%state%v(i,j,1,k1,n0) ! contra v2 = elem(ie)%state%v(i,j,2,k1,n0) ! contra vlatlon(i,j,1) =elem(ie)%D(i,j,1,1)*v1 + elem(ie)%D(i,j,1,2)*v2 ! contra->latlon vlatlon(i,j,2) =elem(ie)%D(i,j,2,1)*v1 + elem(ie)%D(i,j,2,2)*v2 ! contra->latlon E(i,j) = 0.5D0*(vlatlon(i,j,1)**2 + vlatlon(i,j,2)**2) hstar = (elem(ie)%state%p(i,j,1,n0) + pmean)/g hs = elem(ie)%state%ps(i,j)/g h = hstar + hs kenergy(i,j,ie) = hstar*E(i,j) penergy(i,j,ie) = 0.5D0*g*(h**2 - hs**2) diss_p(i,j,k1,ie) = nu_s*(E(i,j)+g*h)*diss_p(i,j,k1,ie)/g + & nu*hstar*(vlatlon(i,j,1)*diss_v(i,j,1,k1,ie)+ & vlatlon(i,j,2)*diss_v(i,j,2,k1,ie)) ! weak dissipation operator has mass matrix "built in" ! integral will multiply by mass matrix again, so lets remove it: ! note: this is not the same as multiplying by rspheremp()! diss_p(i,j,k1,ie)=diss_p(i,j,k1,ie)/elem(ie)%spheremp(i,j) !hs etc. is used only from the first level !but mass conservation output will be done for all levels do k=1,nlev mass(i,j,k,ie) = (elem(ie)%state%p(i,j,k,n0) + pmean)/g enddo end do end do ! global integral < f*g > is computed correctly if either f or g is continious ! enstropy = < PV**2 >, so we need to make the PV continious before computing global integral div(:,:,ie) = divergence_sphere(vlatlon,deriv,elem(ie)) ! latlon vector -> scalar vor(:,:,ie) = vorticity_sphere(vlatlon,deriv,elem(ie)) ! latlon vector -> scalar penst(:,:,ie) = vor(:,:,ie)*elem(ie)%spheremp(:,:) kptr=0 call edgeVpack(edge3, penst(1,1,ie), 1, kptr,ie) end do call bndry_exchangeV(hybrid,edge3) do ie=nets,nete kptr=0 call edgeVunpack(edge3, penst(1,1,ie), 1, kptr, ie) penst(:,:,ie) = penst(:,:,ie)*elem(ie)%rspheremp(:,:) do j=1,np do i=1,np pv(i,j,ie) = penst(i,j,ie) + elem(ie)%fcor(i,j) hstar = (elem(ie)%state%p(i,j,1,n0) + pmean)/g ! For swtc1 p = -pmean and penst is not well defined if (compute_penst) then penst(i,j,ie) = (0.5D0*(penst(i,j,ie) + elem(ie)%fcor(i,j))**2) / hstar else penst(i,j,ie) = 0.0D0 endif end do end do end do do k=1,nlev Imass(k) = global_integral(elem,mass(:,:,k,nets:nete),hybrid,np,nets,nete) enddo Ikenergy = global_integral(elem,kenergy(:,:,nets:nete),hybrid,np,nets,nete) Ipenergy = global_integral(elem,penergy(:,:,nets:nete),hybrid,np,nets,nete) Ienergy = Ikenergy + Ipenergy Ipenst = global_integral(elem,penst(:,:,nets:nete),hybrid,np,nets,nete) Ipv = global_integral(elem,pv(:,:,nets:nete),hybrid,np,nets,nete) Idiv = global_integral(elem,div(:,:,nets:nete),hybrid,np,nets,nete) dissE = global_integral(elem,diss_p(:,:,1,nets:nete),hybrid,np,nets,nete) if (hypervis_order==2) dissE=-dissE ! biharmonic is applied with a minus sign rms_div = global_integral(elem,div(:,:,nets:nete)**2,hybrid,np,nets,nete) rms_vor = global_integral(elem,vor(:,:,nets:nete)**2,hybrid,np,nets,nete) ! planatery vorticity for normalization do ie=nets,nete vor(:,:,ie)=elem(ie)%fcor(:,:) enddo rms_f = global_integral(elem,vor(:,:,nets:nete)**2,hybrid,np,nets,nete) time = Time_at(tl%nstep) if (hybrid%masterthread) then if (tl%nstep==0) then Imass_init = Imass Ienergy_init = Ienergy Ipenst_init = Ipenst Ipv_init = Ipv Idiv_init = Idiv end if write(iounit+0,*)time/secpday,Imass(1) write(iounit+1,*)time/secpday,Ikenergy write(iounit+2,*)time/secpday,Ipenergy write(iounit+3,*)time/secpday,Ipenst write(iounit+4,*)time/secpday,Ipv write(iounit+5,*)time/secpday,Idiv write(iounit+6,*)time/secpday,sqrt(rms_div/rms_f) write(iounit+7,*)time/secpday,sqrt(rms_vor/rms_f) !write(6,*)time/secpday,"mass =",(Imass-Imass_init)/Imass_init !write(6,*)time/secpday,"total energy =",(Ienergy-Ienergy_init)/Ienergy_init !write(6,*)time/secpday,"pot enstrophy =",(Ipenst-Ipenst_init)/Ipenst_init if (tl%nstep>0) then dKE = (Ikenergy-Ikenergy_last)/(time-time_last)/(Ikenergy+Ipenergy) dPE = (Ipenergy-Ipenergy_last)/(time-time_last)/(Ikenergy+Ipenergy) dENS = (Ipenst-Ipenst_last)/(time-time_last)/Ipenst do k=1,nlev dm = (Imass(k)-Imass_last(k))/(time-time_last)/Imass(k) write(6,'(a,i3,a,e11.4,a,e13.6,a,e13.6)') "level =",k," M-M0/M0 =",(Imass(k)-Imass_init(k))/Imass_init(k), & " dM/dt /M = ",dm enddo if (test_case(1:5) /= "swtc1") then write(6,'(a,e11.4,a,e13.6,a,e13.6)') "ENS-ENS0/ENS0=",(Ipenst-Ipenst_init)/Ipenst_init, " dENS/dt /ENS= ",dENS write(6,'(a,e11.4,a,e13.6,a,e11.4,a,e11.4,a)') "E-E0/E0 =",(Ienergy-Ienergy_init)/Ienergy_init,& " dE/dt /E = ",dKE+dPE," (",dKE,",",dPE,")" write(6,'(a,e13.6,f10.3)') " Dissipation/E = ",dissE/(Ikenergy+Ipenergy) endif dm = (Ipv-Ipv_last)/(time-time_last) write(6,'(a,e11.4,a,e13.6,a,e13.6,a,e13.6)') "PV-PV0 =",(Ipv-Ipv_init),& " dPV/dt = ",dm," |vor|/|f|=",sqrt(rms_vor/rms_f) dm = (Idiv-Idiv_last)/(time-time_last) write(6,'(a,e11.4,a,e13.6,a,e13.6,a,e13.6)') "DV-DV0 =",(Idiv-Idiv_init),& " dDV/dt = ",dm," |div|/|f|=",sqrt(rms_div/rms_f) endif end if Imass_last=Imass Ipenergy_last=Ipenergy Ikenergy_last=Ikenergy Ipenst_last=Ipenst Ipv_last=Ipv Idiv_last=Idiv time_last = time end subroutine sweq_invariants ! =========================================== ! tc1_init_pmean ! =========================================== function tc1_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean pmean = gh0_tc1 ! pmean is needed for semi-implicit ! for explicit, it does not matter and it gives negative values ! with monotone limiter of order -1e-11 for some reason if (integration == "explicit" .and. limiter_option /= 0) pmean = 0 if (integration == "runge_kutta" .and. limiter_option /= 0) pmean = 0 if (integration == "full_imp" .and. limiter_option /= 0) pmean = 0 end function tc1_init_pmean ! =========================================== ! tc2_init_pmean ! =========================================== function tc2_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean pmean = gh0_tc2 end function tc2_init_pmean ! =========================================== ! tc5_init_pmean ! =========================================== function tc5_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean pmean = gh0_tc5 end function tc5_init_pmean ! =========================================== ! tc6_init_pmean ! =========================================== function tc6_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean ! gh0_tc6 = g*h0_tc6 h0_tc6=8d3 pmean = gh0_tc6 if (sub_case > 1 ) pmean = g*h0_tc6*sub_case if (sub_case < 0 ) then R_rh = 32 pmean = g*h0_tc6*abs(sub_case) end if end function tc6_init_pmean ! =========================================== ! tc8_init_pmean ! =========================================== function tc8_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean pmean = gh0_tc8 end function tc8_init_pmean ! =========================================== ! sj1_init_pmean ! =========================================== function sj1_init_pmean() result(pmean) implicit none real (kind=real_kind) :: pmean pmean = gh0_sj1 end function sj1_init_pmean ! =========================================== ! ! tc1_init_state: ! ! Initialize for shallow water testcase 1 ! ! =========================================== subroutine tc1_init_state(elem, nets,nete,pmean) type(element_t), intent(inout) :: elem(:) integer, intent(in) :: nets integer, intent(in) :: nete real (kind=real_kind) :: pmean ! ------ where is intent(in) for pmean??? OG ! Local variables integer ie,k,case integer :: nm1 integer :: n0 integer :: np1 ! For swtc1 p = -pmean and penst is not well defined compute_penst = .false. nm1= 1 n0 = 2 np1= 3 pmean = tc1_init_pmean() ! set pmean do ie=nets,nete do k=1,nlev case = sub_case if (sub_case==1 .and. k>1) case=k ! run one sub_case in each level elem(ie)%state%p(:,:,k,n0)=tc1_geopotential(elem(ie)%spherep(:,:),pmean,latc0,lonc0,case) ! elem(ie)%state%p(:,:,k,n0)=0.0D0 elem(ie)%state%p(:,:,k,nm1)=elem(ie)%state%p(:,:,k,n0) elem(ie)%state%p(:,:,k,np1)=0.0D0 elem(ie)%state%v(:,:,:,k,n0)=tc1_velocity(elem(ie)%spherep(:,:),elem(ie)%Dinv) elem(ie)%state%v(:,:,:,k,nm1)=elem(ie)%state%v(:,:,:,k,n0) elem(ie)%state%v(:,:,:,k,np1)=0.0D0 end do elem(ie)%state%ps(:,:)=0.0D0 elem(ie)%state%gradps(:,:,:)=0.0D0 end do end subroutine tc1_init_state ! =========================================== ! ! tc1_velocity: ! ! Reset velocities to initial values ! for a single layer of one element for ! test case 1, during the timestep. ! ! =========================================== function tc1_velocity(sphere,D) result(v) type (spherical_polar_t), intent(in) :: sphere(np,np) real (kind=real_kind), intent(in) :: D(np,np,2,2) real (kind=real_kind) :: v(np,np,2) ! Local variables real (kind=real_kind) :: snlon(np,np) real (kind=real_kind) :: cslon(np,np) real (kind=real_kind) :: snlat(np,np) real (kind=real_kind) :: cslat(np,np) real (kind=real_kind) :: csalpha real (kind=real_kind) :: snalpha real (kind=real_kind) :: V1,V2 integer i,j,k csalpha = COS(alpha) snalpha = SIN(alpha) do j=1,np do i=1,np snlat(i,j) = SIN(sphere(i,j)%lat) cslat(i,j) = COS(sphere(i,j)%lat) snlon(i,j) = SIN(sphere(i,j)%lon) cslon(i,j) = COS(sphere(i,j)%lon) end do end do do k=1,nlev do j=1,np do i=1,np V1 = u0*(cslat(i,j)*csalpha + snlat(i,j)*cslon(i,j)*snalpha) V2 = -u0*(snlon(i,j)*snalpha) ! ===================================================== ! map sphere velocities onto the contravariant cube velocities ! using the D^-T mapping matrix (see Loft notes for details) ! ===================================================== v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end do end function tc1_velocity ! !****************************************************************************** ! ! Subroutines for rotation in most general settings, from 3d code ! !****************************************************************************** ! SUBROUTINE regrot(pxreg,pyreg,pxrot,pyrot,pxcen,pycen,kcall) IMPLICIT NONE ! !---------------------------------------------------------------------- ! !* conversion between regular and rotated spherical coordinates. !* !* pxreg longitudes of the regular coordinates !* pyreg latitudes of the regular coordinates !* pxrot longitudes of the rotated coordinates !* pyrot latitudes of the rotated coordinates !* all coordinates given in degrees n (negative for s) !* and degrees e (negative values for w) !* pxcen regular longitude of the south pole of the rotated grid !* pycen regular latitude of the south pole of the rotated grid !* !* kcall=-1: find regular as functions of rotated coordinates. !* kcall= 1: find rotated as functions of regular coordinates. ! !----------------------------------------------------------------------- ! integer kxdim,kydim,kx,ky,kcall real(kind=real_kind) :: pxreg,pyreg,& pxrot,pyrot,& pxcen,pycen ! !----------------------------------------------------------------------- ! real(kind=real_kind) zsycen,zcycen,zxmxc,zsxmxc,zcxmxc,zsyreg,zcyreg, & zsyrot,zcyrot,zcxrot,zsxrot,zpi,zpih integer jy,jx zpih = DD_PI*0.5d0 ! !---------------------------------------------------------------------- ! zsycen = SIN((pycen+zpih)) zcycen = COS((pycen+zpih)) ! IF (kcall.eq.1) then ! zxmxc = pxreg - pxcen zsxmxc = SIN(zxmxc) zcxmxc = COS(zxmxc) zsyreg = SIN(pyreg) zcyreg = COS(pyreg) zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc zsyrot = max(zsyrot,-1.0D0) zsyrot = min(zsyrot,+1.0D0) ! pyrot = ASIN(zsyrot) ! zcyrot = COS(pyrot) zcxrot = (zcycen*zcyreg*zcxmxc +zsycen*zsyreg)/zcyrot zcxrot = max(zcxrot,-1.0D0) zcxrot = min(zcxrot,+1.0D0) zsxrot = zcyreg*zsxmxc/zcyrot ! pxrot = ACOS(zcxrot) ! IF (zsxrot<0.0) pxrot = -pxrot ! ELSEIF (kcall.eq.-1) then ! zsxrot = SIN(pxrot) zcxrot = COS(pxrot) zsyrot = SIN(pyrot) zcyrot = COS(pyrot) zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot zsyreg = max(zsyreg,-1.0D0) zsyreg = min(zsyreg,+1.0D0) ! pyreg = ASIN(zsyreg) ! zcyreg = COS(pyreg) zcxmxc = (zcycen*zcyrot*zcxrot -& zsycen*zsyrot)/zcyreg zcxmxc = max(zcxmxc,-1.0D0) zcxmxc = min(zcxmxc,+1.0D0) zsxmxc = zcyrot*zsxrot/zcyreg zxmxc = ACOS(zcxmxc) IF (zsxmxc<0.0) zxmxc = -zxmxc ! pxreg = zxmxc + pxcen ! ELSE WRITE(6,'(1x,''invalid kcall in regrot'')') STOP ENDIF END SUBROUTINE regrot SUBROUTINE turnwi(puarg,pvarg,pures,pvres, & pxreg,pyreg,pxrot,pyrot, & pxcen,pycen,kcall) ! IMPLICIT NONE ! !----------------------------------------------------------------------- ! !* turn horizontal velocity components between regular and !* rotated spherical coordinates. ! !* puarg : input u components !* pvarg : input v components !* pures : output u components !* pvres : output v components !* pa : transformation coefficients !* pb : -"- !* pc : -"- !* pd : -"- !* pxreg : regular longitudes !* pyreg : regular latitudes !* pxrot : rotated longitudes !* pyrot : rotated latitudes !* kxdim : dimension in the x (longitude) direction !* kydim : dimension in the y (latitude) direction !* kx : number of gridpoints in the x direction !* ky : number of gridpoints in the y direction !* pxcen : regular longitude of the south pole of the !* transformed grid !* pycen : regular latitude of the south pole of the !* transformed grid !* !* kcall < 0 : find wind components in regular coordinates !* from wind components in rotated coordinates !* kcall > 0 : find wind components in rotated coordinates !* from wind components in regular coordinates !* note that all coordinates are given in degrees n and degrees e. !* (negative values for s and w) ! !----------------------------------------------------------------------- integer kxdim,kydim,kx,ky,kcall real(kind=real_kind) puarg,pvarg, & pures,pvres, & pa, pb, & pc, pd, & pxreg,pyreg, & pxrot,pyrot real(kind=real_kind) pxcen,pycen ! !----------------------------------------------------------------------- ! integer jy,jx real(kind=real_kind) zpih,zsyc,zcyc,zsxreg,zcxreg,zsyreg,zcyreg,zxmxc,& zsxmxc,zcxmxc,zsxrot,zcxrot,zsyrot,zcyrot ! !----------------------------------------------------------------------- ! IF (kcall.eq.1) then zpih = DD_PI*0.5d0 zsyc = SIN(pycen+zpih) zcyc = COS(pycen+zpih) ! zsxreg = SIN(pxreg) zcxreg = COS(pxreg) zsyreg = SIN(pyreg) zcyreg = COS(pyreg) ! zxmxc = pxreg - pxcen zsxmxc = SIN(zxmxc) zcxmxc = COS(zxmxc) ! zsxrot = SIN(pxrot) zcxrot = COS(pxrot) zsyrot = SIN(pyrot) zcyrot = COS(pyrot) ! pa = zcyc*zsxmxc*zsxrot + zcxmxc*zcxrot pb = zcyc*zcxmxc*zsyreg*zsxrot - zsyc*zcyreg*zsxrot - & zsxmxc*zsyreg*zcxrot pc = zsyc*zsxmxc/zcyrot pd = (zsyc*zcxmxc*zsyreg + zcyc*zcyreg)/zcyrot ! pures = pa*puarg + pb*pvarg pvres = pc*puarg + pd*pvarg ELSEIF (kcall.eq.-1) then zpih = DD_PI*0.5d0 zsyc = SIN(pycen+zpih) zcyc = COS(pycen+zpih) ! zsxreg = SIN(pxreg) zcxreg = COS(pxreg) zsyreg = SIN(pyreg) zcyreg = COS(pyreg) ! zxmxc = pxreg - pxcen zsxmxc = SIN(zxmxc) zcxmxc = COS(zxmxc) ! zsxrot = SIN(pxrot) zcxrot = COS(pxrot) zsyrot = SIN(pyrot) zcyrot = COS(pyrot) ! pa = zcxmxc*zcxrot + zcyc*zsxmxc*zsxrot pb = zcyc*zsxmxc*zcxrot*zsyrot + zsyc*zsxmxc*zcyrot -& zcxmxc*zsxrot*zsyrot pc =-zsyc*zsxrot/zcyreg pd = (zcyc*zcyrot - zsyc*zcxrot*zsyrot)/zcyreg ! pures = pa*puarg + pb*pvarg pvres = pc*puarg + pd*pvarg ELSE write(6,'(1x,''invalid kcall in turnwi'')') STOP ENDIF END SUBROUTINE turnwi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! =========================================== ! ! tc1_geopotential: ! ! Initialize cosine bell centered at latc, lonc ! ! =========================================== function tc1_geopotential(sphere, pmean, latc, lonc,k) result(p) type (spherical_polar_t), intent(in) :: sphere(np,np) ! spherical coordinates of element real (kind=real_kind), intent(in) :: pmean ! mean geopotential real (kind=real_kind), intent(in) :: latc ! latitude of center of cosine bell real (kind=real_kind), intent(in) :: lonc ! longitude of center of cosine bell real (kind=real_kind) :: p(np,np) ! geopotential real (kind=real_kind) :: snlat,cslat,lon,lat real (kind=real_kind) :: h real (kind=real_kind) :: r,A real (kind=real_kind) :: h2,lon_min,lon_max integer i,j,k h2 = 0.25D0*rr lon_min = lonc - h2 lon_max = lonc + h2 do j=1,np do i=1,np snlat=SIN(sphere(i,j)%lat) cslat=COS(sphere(i,j)%lat) lon =sphere(i,j)%lon lat =sphere(i,j)%lat A =SIN(latc)*snlat + COS(latc)*cslat*COS(lon - lonc) ! face 4 this == - cslat*snlon r=ACOS(A) if (r lon_min).and.(lon -DD_PI).and.(lat 2.0D0*DD_PI + latc - rr))then p(i,j)=0.0D0 end if end if p(i,j)=g*p(i,j) - pmean end if end do end do end function tc1_geopotential ! =========================================== ! ! tc1_phi: (nair) for unstaggered grids ! ! Initialize cosine bell centered at latc, lonc on (np,np) points ! ! =========================================== function tc1_phi(sphere, pmean, latc, lonc) result(ppsi) type (spherical_polar_t), intent(in) :: sphere(np,np) ! spherical coordinates of element real (kind=real_kind), intent(in) :: pmean ! mean geopotential real (kind=real_kind), intent(in) :: latc ! latitude of center of cosine bell real (kind=real_kind), intent(in) :: lonc ! longitude of center of cosine bell real (kind=real_kind) :: ppsi(np,np) ! geopotential real (kind=real_kind) :: snlat,cslat,lon, cslon, coef real (kind=real_kind) :: h, gh00, trm real (kind=real_kind) :: r,A, csalpha, snalpha integer i,j csalpha = COS(alpha) snalpha = SIN(alpha) do j=1,np do i=1,np snlat=SIN(sphere(i,j)%lat) cslat=COS(sphere(i,j)%lat) cslon=COS(sphere(i,j)%lon) lon =sphere(i,j)%lon A = SIN(latc)*snlat + COS(latc)*cslat*COS(lon - lonc) ! face 4 this == - cslat*snlon r=ACOS(A) if (r=6) toy_chemistry=1 endif ! original units: bad mix of dimensionless and dimensional: ! Tperiod=5*24*60*60 !for 5 days ! dimensional unts. Run 12 days: Tperiod=12*24*60*60 !for 12 days Kcoef=Kcoef*rearth*5/Tperiod nm1= 1 n0 = 2 np1= 3 do ie=nets,nete do k=1,nlev elem(ie)%state%p(:,:,k,n0)=swirl_init_tracer(elem(ie)%spherep(:,:),k) elem(ie)%state%v(:,:,:,k,n0)=swirl_velocity(0.0D0,elem(ie)%spherep(:,:),elem(ie)%Dinv) elem(ie)%state%v(:,:,:,k,nm1)=elem(ie)%state%v(:,:,:,k,n0) elem(ie)%state%v(:,:,:,k,np1)=elem(ie)%state%v(:,:,:,k,n0) elem(ie)%state%p(:,:,k,nm1)=elem(ie)%state%p(:,:,k,n0) elem(ie)%state%p(:,:,k,np1)=elem(ie)%state%p(:,:,k,n0) end do elem(ie)%state%ps(:,:)=0.0D0 elem(ie)%state%gradps(:,:,:)=0.0D0 end do g_sw_output=1.0D0 if((kmass_swirl<0).or.(kmass_swirl>nlev))then else kmass=kmass_swirl endif end subroutine swirl_init_state ! =========================================== ! ! swirl_init_tracer: ! ! Initializes many tracers ! ! =========================================== recursive function swirl_init_tracer(sphere,level) result(p) type (spherical_polar_t), intent(in) :: sphere(np,np) ! spherical coordinates of element integer, intent(in) :: level real (kind=real_kind) :: p(np,np) ! geopotential real (kind=real_kind) :: snlat,cslat,lon,lat, lon_min, lon_max real (kind=real_kind) :: h1,h2 real (kind=real_kind) :: r1,r2,A1,A2, bellrad, latc, lonc,xc,yc,zc,snlon,cslon,x1,y1,z1 real (kind=real_kind) :: lon_1,lat_1,lon_2,lat_2 real (kind=real_kind) :: k1,k2,r,D,Cly,lat_center,lon_center,density,dminusr(np,np) integer i,j,k ! ! Nick Lopez 7/2015 added terminator toy chemistry ! ! layers in swirl: ! lev=1: COS BELLS ! lev=2: GAUSS BELLS ! lev=3: SLOTTED CYLINDERS ! lev=4: CONST=1 ! lev=5: "Cl" from Lauritzen et al. Toy Chemistry test ! lev=6: "Cl2" from Lauritzen et al. Toy Chemistry test ! Note: chemistry is applied to tracers 5 & 6, if they are present !------------------------ levels below are more or less temporary ! lev=7: Tracer 6 + Tracer 2, to test if linear relations are preserved ! lev=8: metric for level 1 ! lev=9: metric for level 2 ! lev=10: metric for level 3 ! lev=11: temporary (right now is used to mark elements) ! lev=12: is filter applied to level 1 and where ! lev=13: is filter applied to level 2 and where ! lev=14: is filter applied to level 3 and where ! code used for tracers 5 and 6 lat_center = 1.0d0/9.0d0*DD_PI lon_center = 2.0d0/3.0d0*DD_PI Cly = 4.0d-6 do i=1, np do j=1, np k1 = SIN(sphere(i,j)%lat)*SIN(lat_center) + COS(sphere(i,j)%lat)*COS(lat_center)*COS(sphere(i,j)%lon-lon_center) if (k1<0.0d0) then k1 = 0.0d0 endif k2 = 1.0d0 r = k1/(4.0d0*k2) D = sqrt(r**2.0d0 + 2.0d0*r*Cly) dminusr(i,j) = ( D-r ) enddo enddo if(level==1)then bellrad=bellradius !bellradius=1/2 !these values for true swirl test cases lon_1=lon1; lat_1=lat1;lon_2=lon2; lat_2=lat2; !these values below make bells apart by more than 1 element if NE=20 !lon_1=3.0*DD_PI/4.0; lat_1=0.0;lon_2=5.0*DD_PI/4.0; lat_2=0.0; do j=1,np do i=1,np snlat=SIN(sphere(i,j)%lat) cslat=COS(sphere(i,j)%lat) lon =sphere(i,j)%lon lat =sphere(i,j)%lat A1=SIN(lat_1)*snlat + COS(lat_1)*cslat*COS(lon - lon_1) A2=SIN(lat_2)*snlat + COS(lat_2)*cslat*COS(lon - lon_2) r1=ACOS(A1) r2=ACOS(A2) if (r1=h2)) then p(i,j)=tracer_highest+tracer_lowest endif if ((r1=h2)) then p(i,j)=tracer_highest+tracer_lowest endif if ((r15.d0/12.d0*bellrad)) then p(i,j)=tracer_highest+tracer_lowest endif enddo enddo !initializing density elseif(level==4)then p(:,:)=1.0D0 elseif(level==5)then p(:,:) = dminusr(:,:) ! Tracer "Cl" elseif(level==6)then p(:,:) = 0.5d0*(Cly-dminusr(:,:)) ! tracer "Cl_2" elseif(level==7)then p(:,:)=swirl_init_tracer(sphere,6)+swirl_init_tracer(sphere,2) !this level is for observing metric of level 1 elseif(level==8)then p(:,:)=0.0d0 !this level is for observing metric of level 2 elseif(level==9)then p(:,:)=0.0d0 !this level is for observing metric of level 3 elseif(level==10)then p(:,:)=0.0d0 !this is another level for debugging or for smoothness measure !or for int smooth cyllinder elseif(level==11)then do j=1,np do i=1,np !lon_1=lon1; lat_1=lat1;lon_2=lon2; lat_2=lat2; lon_1=3.0*DD_PI/4.0; lat_1=0.0;lon_2=5.0*DD_PI/4.0; lat_2=0.0; snlat=SIN(sphere(i,j)%lat) cslat=COS(sphere(i,j)%lat) snlon=SIN(sphere(i,j)%lon) cslon=COS(sphere(i,j)%lon) xc = cos(lat_2)*cos(lon_2) yc = cos(lat_2)*sin(lon_2) zc = sin(lat_2) x1 = cslat*cslon y1 = cslat*snlon z1 = snlat r1 = (x1-xc)**2 + (y1-yc)**2 + (z1-zc)**2 h1 = 0.5d0+0.5d0*tanh(8.0d0*r1+3.d0)*tanh(3.d0-8.0d0*r1) p(i,j)=h1 end do end do !this level is for observing limiter applied to level 1 elseif(level==12)then p(:,:)=0.0d0 !this level is for observing limiter applied to level 2 elseif(level==13)then p(:,:)=0.0d0 !this level is for observing limiter applied to level 3 elseif(level==14)then p(:,:)=0.0d0 endif end function swirl_init_tracer !-------------------------------------------------------------------- ! velocities for test SWIRL ! ! modif by OG !-------------------------------------------------------------------- function swirl_velocity(t,sphere,D) result(v) type (spherical_polar_t), intent(in) :: sphere(np,np) real (kind=real_kind), intent(in) :: t,D(np,np,2,2) real (kind=real_kind) :: v(np,np,2) ! Local variables real (kind=real_kind) :: snhalflon(np,np) real (kind=real_kind) :: snlon(np,np) real (kind=real_kind) :: sndoublelat(np,np) real (kind=real_kind) :: sndoublelon(np,np) real (kind=real_kind) :: cslat(np,np) real (kind=real_kind) :: csperiod, PI real (kind=real_kind) :: V1,V2 integer i,j,k real (kind=real_kind) :: sphlon,sphlat real (kind=real_kind) :: v1_old,v2_old,v1_new,v2_new,lonnew,latnew,shift_for_rotation real (kind=real_kind) :: vel_pure_rotation_marker if(add_pure_rotation)then shift_for_rotation=t*DD_PI/Tperiod vel_pure_rotation_marker=1.0D0 else shift_for_rotation=0.0D0 vel_pure_rotation_marker=0.0D0 endif csperiod=COS(t*DD_PI/Tperiod) if(sub_case==1)then do j=1,np do i=1,np !peter's code for adding zonal flow ! tt = 5 ! ck = 2 ! slon = sin(lon-(t)*pi2/tt)!solid-body rotation added ! slon2 = sin(two*(lon-(t)*pi2/tt))!solid-body rotation added ! clon = cos(lon-(t)*pi2/tt)!solid-body rotation added ! clon2 = cos(two*(lon-(t)*pi2/tt))!solid-body rotation added ! sslm = (sin(half*(lon-(t)*pi2/tt)))**2 ! ! uexact = ck*slon*slon*SIN(two*lat)*COS(t*omega) + clat*pi2/(tt) ! vexact = ck*slon2*clat*COS(t*omega) ! udc = two*ck*slon*slon*slat*COS(t*omega) + pi2/tt sphlon=sphere(i,j)%lon-shift_for_rotation*2.0d0 sphlat=sphere(i,j)%lat CALL regrot(sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,1) snlon(i,j) = SIN(lonnew) snhalflon(i,j) = SIN(lonnew/2.0D0) cslat(i,j) = COS(latnew) sndoublelat(i,j) = SIN(latnew*2.0D0) sndoublelon(i,j) = SIN(lonnew*2.0D0) do k=1,nlev V1 = Kcoef*(snhalflon(i,j))**2*sndoublelat(i,j)*csperiod & +cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker V2 = Kcoef/2.0D0*snlon(i,j)*cslat(i,j)*csperiod v1_old=V1 v2_old=V2 CALL turnwi(v1_old,v2_old,& v1_new,v2_new,sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,-1) IF (ABS(DD_PI*0.5d0-sphlat)<1.0E-8.OR.ABS(DD_PI*0.5d0+sphlat)<1.0E-8) THEN v2_new = 0.0d0 endif IF (ABS(v1_new)<1.0E-10) v1_new=0.0d0 V1=v1_new V2=v2_new v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end do elseif(sub_case==2)then do j=1,np do i=1,np sphlon=sphere(i,j)%lon-shift_for_rotation*2 sphlat=sphere(i,j)%lat CALL regrot(sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,1) snlon(i,j) = SIN(lonnew) snhalflon(i,j) = SIN(lonnew/2.0D0) cslat(i,j) = COS(latnew) sndoublelat(i,j) = SIN(latnew*2.0D0) sndoublelon(i,j) = SIN(lonnew*2.0D0) do k=1,nlev V1 = Kcoef*(snlon(i,j))**2*sndoublelat(i,j)*csperiod +& cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker V2 = Kcoef*sndoublelon(i,j)*cslat(i,j)*csperiod v1_old=V1 v2_old=V2 CALL turnwi(v1_old,v2_old,& v1_new,v2_new,sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,-1) IF (ABS(DD_PI*0.5d0-sphlat)<1.0E-8.OR.ABS(DD_PI*0.5d0+sphlat)<1.0E-8) THEN v2_new = 0.0d0 endif IF (ABS(v1_new)<1.0E-10) v1_new=0.0d0 V1=v1_new V2=v2_new v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end do elseif(sub_case==3)then do j=1,np do i=1,np sphlon=sphere(i,j)%lon-shift_for_rotation*2.0d0 sphlat=sphere(i,j)%lat CALL regrot(sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,1) snlon(i,j) = SIN(lonnew) snhalflon(i,j) = SIN(lonnew/2.0D0) cslat(i,j) = COS(latnew) sndoublelat(i,j) = SIN(latnew*2.0D0) sndoublelon(i,j) = SIN(lonnew*2.0D0) do k=1,nlev V1 = -Kcoef*(snhalflon(i,j))**2*sndoublelat(i,j)*cslat(i,j)**2*csperiod & +cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker V2 = Kcoef/2.0D0*snlon(i,j)*cslat(i,j)**3*csperiod v1_old=V1 v2_old=V2 CALL turnwi(v1_old,v2_old,& v1_new,v2_new,sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,-1) IF (ABS(DD_PI*0.5d0-sphlat)<1.0E-8.OR.ABS(DD_PI*0.5d0+sphlat)<1.0E-8) THEN v2_new = 0.0d0 endif IF (ABS(v1_new)<1.0E-10) v1_new=0.0d0 V1=v1_new V2=v2_new v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end do elseif(sub_case==4)then do j=1,np do i=1,np sphlon=sphere(i,j)%lon-shift_for_rotation*2.0d0 sphlat=sphere(i,j)%lat CALL regrot(sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,1) snlon(i,j) = SIN(lonnew) snhalflon(i,j) = SIN(lonnew/2.0D0) cslat(i,j) = COS(latnew) sndoublelat(i,j) = SIN(latnew*2.0D0) sndoublelon(i,j) = SIN(lonnew*2.0D0) do k=1,nlev !true case 4 V1 = Kcoef*(snlon(i,j))**2*sndoublelat(i,j)*csperiod & +cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker V2 = Kcoef*sndoublelon(i,j)*cslat(i,j)*csperiod !experimental case 4 with flow 'away' from poles ! V1 = 0.5*Kcoef*(snlon(i,j))**2*sndoublelat(i,j)*csperiod & ! +cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker ! V2 = 0.5*Kcoef*sndoublelon(i,j)*cslat(i,j)*csperiod !only rotation to test limiters ! V1 = cslat(i,j)*2.0D0*DD_PI*rearth/Tperiod*vel_pure_rotation_marker ! V2 = 0.0d0 v1_old=V1 v2_old=V2 CALL turnwi(v1_old,v2_old,& v1_new,v2_new,sphlon,sphlat,lonnew,latnew,0.0d0,-0.5d0*DD_PI+rotangle,-1) IF (ABS(DD_PI*0.5d0-sphlat)<1.0E-8.OR.ABS(DD_PI*0.5d0+sphlat)<1.0E-8) THEN v2_new = 0.0d0 endif IF (ABS(v1_new)<1.0E-10) v1_new=0.0d0 V1=v1_new V2=v2_new v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end do endif end function swirl_velocity ! =================================================================== ! ! sj1_init_state: ! ! Initialize strong jet case 1: Galewski et al. ! ! =================================================================== subroutine sj1_init_state(elem,nets,nete,hybrid,pmean,deriv) type(element_t), target,intent(inout) :: elem(:) integer, intent(in) :: nets integer, intent(in) :: nete type (hybrid_t) :: hybrid real (kind=real_kind) :: pmean type (derivative_t) :: deriv ! Local variables real (kind=real_kind) :: a,b real (kind=real_kind) :: mean_balance real (kind=real_kind) :: pert(np,np) real (kind=real_kind) :: phi_balance(np,np,nets:nete) type (quadrature_t) :: gs integer :: ie,k integer :: nm1 integer :: n0 integer :: np1 type (element_t), pointer :: pElem,sElem logical :: Debug=.FALSE. nm1= 1 n0 = 2 np1= 3 a =-DD_PI/2.0D0 b = DD_PI/2.0D0 gs=gauss(Ngs_sj1) Ibalance_tot = gaussian_int(sj1_integrand,a,b,gs) print *,"Ibalance_tot=",Ibalance_tot deallocate(gs%points) deallocate(gs%weights) print *,"done..." do ie=nets,nete pElem => elem(ie) phi_balance(:,:,ie) = sj1_balance(pElem%spherep(:,:),np) end do pmean = g*h0_sj1 gh0_sj1 = pmean mean_balance = global_integral(elem,phi_balance(:,:,nets:nete),hybrid,np,nets,nete) if(Debug)print *,"BEFORE mean height=",mean_balance/g if(Debug)print *,"BEFORE pmean=",pmean/g pmean = gh0_sj1 + mean_balance gh0_sj1 = pmean mean_balance = pmean-global_integral(elem,phi_balance(:,:,nets:nete),hybrid,np,nets,nete) if(Debug)print *,"AFTER mean height=",mean_balance/g if(Debug)print *,"AFTER pmean=",pmean/g do ie=nets,nete pElem => elem(ie) pert(:,:) = sj1_perturbation(pElem%spherep(:,:),np) pElem%state%ps(:,:) = 0.0D0 do k=1,nlev pElem%state%p(:,:,k,n0) =pmean - phi_balance(:,:,ie)+ g*pert(:,:) pElem%state%p(:,:,k,nm1)=pElem%state%p(:,:,k,n0) pElem%state%p(:,:,k,np1)=pElem%state%p(:,:,k,n0) pElem%state%v(:,:,:,k,n0)=sj1_velocity_cubedsphere(pElem%spherep(:,:),pElem%Dinv) pElem%state%v(:,:,:,k,nm1)=pElem%state%v(:,:,:,k,n0) pElem%state%v(:,:,:,k,np1)=pElem%state%v(:,:,:,k,n0) end do pElem%state%gradps(:,:,:) = 0.0D0 end do pmean=0.0D0 end subroutine sj1_init_state function sj1_velocity(lat) result(ulat) real (kind=real_kind),intent(in) :: lat real (kind=real_kind) :: ulat,en en = exp(-((2.0D0/(lat1_sj1 - lat0_sj1))**2.0D0)) if((lat0_sj1.lt.lat).AND.(lat.lt.lat1_sj1))then ulat = (u0_sj1/en)*exp(1.0D0/((lat-lat0_sj1)*(lat-lat1_sj1))) else ulat = 0.0_real_kind endif end function sj1_velocity ! ======================================================= ! ! sj1_integrand: ! ! Scaled integrand function needed to integrate the ! geostrophic balance equation using the trapezoid rule. ! Integrand in units of u0_sj1*f*rearth (m^2/s^2) ! ! ======================================================= function sj1_integrand(lat) result(Ival) real (kind=real_kind), intent(in) :: lat real (kind=real_kind) :: Ival real (kind=real_kind) :: ulat ! dimensionless velocity ulat = sj1_velocity(lat) Ival = rearth*ulat*(2.0D0*omega*sin(lat) + tan(lat)*ulat/rearth) end function sj1_integrand ! ========================================================== ! ! sj1_balance: ! ! Compute geostrophic integral balance equation for the ! Barotropic Instability Test Case of Polvani. ! ! =========================================================== function sj1_balance(sphere,npts) result(p) integer, intent(in) :: npts type (spherical_polar_t), intent(in) :: sphere(npts,npts) real (kind=real_kind) :: p(npts,npts) ! Local variables real (kind=real_kind) :: lat real (kind=real_kind) :: lon real (kind=real_kind) :: a, b real (kind=real_kind) :: Ibalance real (kind=real_kind) :: f real (kind=real_kind) :: pi2 ! pi/2 integer i,j type (quadrature_t) :: gs gs=gauss(Ngs_sj1) do j=1,npts do i=1,npts lat = sphere(i,j)%lat lon = sphere(i,j)%lon pi2 = 0.5D0*DD_PI b = lat a = -pi2 Ibalance = gaussian_int(sj1_integrand,a,b,gs) p(i,j) = Ibalance end do end do deallocate(gs%weights) deallocate(gs%points) end function sj1_balance ! =========================================== ! ! sj1_perturbation: ! ! Compute the geopotential perturbation ! from geostrophic balance for the Barotropic ! Instability Test Case of Polvani. ! ! =========================================== function sj1_perturbation(sphere,npts) result(p) integer, intent(in) :: npts type (spherical_polar_t), intent(in) :: sphere(npts,npts) real (kind=real_kind) :: p(npts,npts) ! Local variables real (kind=real_kind) :: lat real (kind=real_kind) :: lon real (kind=real_kind) :: alon, blat real (kind=real_kind) :: secha,sechb real (kind=real_kind) :: ghprime integer i,j do j=1,npts do i=1,npts lat = sphere(i,j)%lat lon = sphere(i,j)%lon if(lon.ge.DD_PI)then lon = lon -2.0D0*DD_PI end if p(i,j) = hhat_sj1*cos(lat)*exp(-((lon*alpha_sj1)**2.0D0 + ((lat2_sj1 - lat)*beta_sj1)**2.0D0 )) end do end do end function sj1_perturbation ! =========================================== ! ! sj1_velocity: ! ! Set initial zonal velocity field for barotropic ! instability test case of Polvani. ! ! =========================================== function sj1_velocity_cubedsphere(sphere,D) result(v) type (spherical_polar_t), intent(in) :: sphere(np,np) real (kind=real_kind), intent(in) :: D(np,np,2,2) real (kind=real_kind) :: v(np,np,2) ! Local variables real (kind=real_kind) :: lat real (kind=real_kind) :: V1 integer i,j do j=1,np do i=1,np lat = sphere(i,j)%lat V1 = sj1_velocity(lat) ! Map onto contravariant velocities v(i,j,1)= V1*D(i,j,1,1) v(i,j,2)= V1*D(i,j,2,1) end do end do end function sj1_velocity_cubedsphere subroutine sj1_errors(elem, iounit, tl, pmean, fstub, simday, hybrid, nets, nete, par) type(element_t), intent(in) :: elem(:) integer , intent(in) :: iounit type (TimeLevel_t) , intent(in) :: tl ! model time struct real (kind=real_kind), intent(in) :: pmean ! mean geopotential character(len=*) , intent(in) :: fstub ! file path stub integer , intent(in) :: simday ! current day of simulation type (hybrid_t) , intent(in) :: hybrid integer , intent(in) :: nets integer , intent(in) :: nete type(parallel_t) , intent(in) :: par ! Local variables real (kind=real_kind) :: pt(np,np,nets:nete) real (kind=real_kind) :: p(np,np,nets:nete) real (kind=real_kind) :: vt(np,np,2,nets:nete) real (kind=real_kind) :: v(np,np,2,nets:nete) real (kind=real_kind) :: l1,l2,linf integer ie,npts if ((tl%nstep == 1).or.(tl%nstep == 0)) then if (hybrid%masterthread) then open(iounit+0,file=TRIM(output_prefix)//"sjtc1.l1.errors",form="formatted") open(iounit+1,file=TRIM(output_prefix)//"sjtc1.l2.errors",form="formatted") open(iounit+2,file=TRIM(output_prefix)//"sjtc1.linf.errors",form="formatted") end if end if do ie=nets,nete p(:,:,ie) = elem(ie)%state%p(:,:,1,tl%n0) + pmean v(:,:,:,ie) = elem(ie)%state%v(:,:,:,1,tl%n0) end do end subroutine sj1_errors subroutine toy_chemistry_forcing(elem, nets, nete, tl, dt) ! ---- Nick Lopez 7/2015 type (element_t) , intent(inout), target :: elem(:) integer , intent(in) :: nets integer , intent(in) :: nete type (TimeLevel_t), intent(in) :: tl real (kind=real_kind), intent(in) :: dt real (real_kind) :: lat_center, lon_center, k1, k2 real (real_kind) :: cl, cl2, cly, r, d, e, el, f_cl, f_cl2 real (real_kind) :: lat,lon integer :: cl_index,cl2_index,i,j,ie,n0 integer, parameter :: allowNegChem=1 ! Center of k1. This is essentially a longitude offset of PI from Fig.1 of Lauritzen et. al lat_center = DD_PI/9 lon_center = 2*DD_PI/3 cl_index = 5 cl2_index = 6 do ie=nets, nete do i=1, np do j=1, np lat = elem(ie)%spherep(i,j)%lat lon = elem(ie)%spherep(i,j)%lon cl = elem(ie)%state%p(i,j,cl_index,tl%n0) cl2 = elem(ie)%state%p(i,j,cl2_index,tl%n0) cly = cl + 2.0d0*cl2 k1 = SIN(lat)*SIN(lat_center) + COS(lat)*COS(lat_center)*COS(lon-lon_center) if (k1<0.0) then k1 = 0.0d0 endif k2 = 1.0d0 r = k1 / (4.0d0*k2) d = sqrt( r*r + 2.0d0*r*cly ) e = exp( -4.0d0*k2*d*dt) if (abs(d*k2*dt) .gt. 1.0d-16) then el = (1.0d0-e) / (d*dt) else el = 4.0d0*k2 endif f_cl = -el * (cl-d+r) * (cl+d+r) / (1.0d0 + e + dt*el*(cl+r)) f_cl2 = -f_cl/2.0d0 if (allowNegChem==0) then if (cl + f_cl*dt < 0.0d0) then f_cl = (0.0d0 - cl) / dt f_cl2 = -f_cl/2.0d0 elem(ie)%state%p(i,j,cl_index,tl%n0) = 0.0d0 elem(ie)%state%p(i,j,cl2_index,tl%n0) = cl2 + f_cl2*dt else if (cl2 + f_cl2*dt < 0.0d0) then f_cl2 = (0.0d0 - cl2) / dt f_cl = -f_cl2*2.0d0 elem(ie)%state%p(i,j,cl_index,tl%n0) = cl + f_cl*dt elem(ie)%state%p(i,j,cl2_index,tl%n0) = 0.0d0 else elem(ie)%state%p(i,j,cl_index,tl%n0) = cl + f_cl*dt elem(ie)%state%p(i,j,cl2_index,tl%n0) = cl2 + f_cl2*dt endif else ! Can allow chemistry to force negative values elem(ie)%state%p(i,j,cl_index,tl%n0) = cl + f_cl*dt elem(ie)%state%p(i,j,cl2_index,tl%n0) = cl2 + f_cl2*dt endif enddo enddo enddo end subroutine ! =========================================== ! ! planar_dbl_vrtx_init_state: ! ! Initialize planar double vortex test case ! ! =========================================== subroutine planar_dbl_vrtx_init_state(elem,nets,nete,pmean,deriv) type(element_t), intent(inout) :: elem(:) integer, intent(in) :: nets integer, intent(in) :: nete real (kind=real_kind) :: pmean type (derivative_t) :: deriv ! Local variables integer ie,k integer :: nm1 integer :: n0 integer :: np1 nm1= 1 n0 = 2 np1= 3 pmean = gh0_dbv do ie=nets,nete elem(ie)%fcor(:,:)= f_dbv elem(ie)%state%ps(:,:)=0.0D0 do k=1,nlev elem(ie)%state%p(:,:,k,n0)=planar_dbl_vrtx_geopotential(elem(ie)%spherep(:,:)) elem(ie)%state%p(:,:,k,nm1)=elem(ie)%state%p(:,:,k,n0) elem(ie)%state%p(:,:,k,np1)=0.0D0 elem(ie)%state%v(:,:,:,k,n0)=planar_dbl_vrtx_velocity(elem(ie)%spherep(:,:),elem(ie)%Dinv) elem(ie)%state%v(:,:,:,k,nm1)=elem(ie)%state%v(:,:,:,k,n0) elem(ie)%state%v(:,:,:,k,np1)=0.0D0 end do end do end subroutine planar_dbl_vrtx_init_state ! =========================================== ! ! planar_dbl_vrtx_velocity: ! ! Reset velocities to initial values ! for a single layer of one element for ! double vortex test case, during the timestep. ! ! =========================================== function planar_dbl_vrtx_velocity(sphere,D) result(v) type (spherical_polar_t), intent(in) :: sphere(np,np) real (kind=real_kind), intent(in) :: D(np,np,2,2) real (kind=real_kind) :: v(np,np,2) ! Local variables real (kind=real_kind) :: V1,V2 real (kind=real_kind) :: xprime1,yprime1,xprime2,yprime2 real (kind=real_kind) :: xprimeprime1,xprimeprime2,yprimeprime1,yprimeprime2 real (kind=real_kind) :: xc1, xc2, yc1, yc2 integer i,j xc1 = (0.5D0 - ox_dbv)*Lx xc2 = (0.5D0 + ox_dbv)*Lx yc1 = (0.5D0 - oy_dbv)*Ly yc2 = (0.5D0 + oy_dbv)*Ly do j=1,np do i=1,np xprime1 = 1.0D0/(DD_PI*sigmax_dbv) * SIN(DD_PI/Lx * (sphere(i,j)%lon - xc1)) yprime1 = 1.0D0/(DD_PI*sigmay_dbv) * SIN(DD_PI/Ly * (sphere(i,j)%lat - yc1)) xprime2 = 1.0D0/(DD_PI*sigmax_dbv) * SIN(DD_PI/Lx * (sphere(i,j)%lon - xc2)) yprime2 = 1.0D0/(DD_PI*sigmay_dbv) * SIN(DD_PI/Ly * (sphere(i,j)%lat - yc2)) xprimeprime1 = 1.0D0/(2.0D0*DD_PI*sigmax_dbv) * SIN(2.0D0*DD_PI/Lx * (sphere(i,j)%lon - xc1)) yprimeprime1 = 1.0D0/(2.0D0*DD_PI*sigmay_dbv) * SIN(2.0D0*DD_PI/Ly * (sphere(i,j)%lat - yc1)) xprimeprime2 = 1.0D0/(2.0D0*DD_PI*sigmax_dbv) * SIN(2.0D0*DD_PI/Lx * (sphere(i,j)%lon - xc2)) yprimeprime2 = 1.0D0/(2.0D0*DD_PI*sigmay_dbv) * SIN(2.0D0*DD_PI/Ly * (sphere(i,j)%lat - yc2)) V1 = -g*dh_dbv/(f_dbv * sigmay_dbv*Ly) * (yprimeprime1 * EXP(-0.5D0 * (xprime1 * xprime1 + yprime1 * yprime1)) + yprimeprime2 * EXP(-0.5D0 * (xprime2 * xprime2 + yprime2 * yprime2))) V2 = g*dh_dbv/(f_dbv * sigmax_dbv*Lx) * (xprimeprime1 * EXP(-0.5D0 * (xprime1 * xprime1 + yprime1 * yprime1)) + xprimeprime2 * EXP(-0.5D0 * (xprime2 * xprime2 + yprime2 * yprime2))) ! ===================================================== ! map physical domain velocities onto the contravariant reference domain velocities ! using the D mapping matrix (see Loft notes for details) ! ===================================================== v(i,j,1)= V1*D(i,j,1,1) + V2*D(i,j,1,2) v(i,j,2)= V1*D(i,j,2,1) + V2*D(i,j,2,2) end do end do end function planar_dbl_vrtx_velocity ! =========================================== ! ! planar_dbl_vrtx_geopotential: ! ! Set geopotential to initial values ! ! =========================================== function planar_dbl_vrtx_geopotential(sphere) result(p) type (spherical_polar_t), intent(in) :: sphere(np,np) real (kind=real_kind) :: p(np,np) ! Local variables real (kind=real_kind) :: xprime1,yprime1,xprime2,yprime2 real (kind=real_kind) :: xc1, xc2, yc1, yc2 integer i,j xc1 = (0.5D0 - ox_dbv)*Lx xc2 = (0.5D0 + ox_dbv)*Lx yc1 = (0.5D0 - oy_dbv)*Ly yc2 = (0.5D0 + oy_dbv)*Ly do j=1,np do i=1,np xprime1 = 1.0D0/(DD_PI*sigmax_dbv) * SIN(DD_PI/Lx * (sphere(i,j)%lon - xc1)) yprime1 = 1.0D0/(DD_PI*sigmay_dbv) * SIN(DD_PI/Ly * (sphere(i,j)%lat - yc1)) xprime2 = 1.0D0/(DD_PI*sigmax_dbv) * SIN(DD_PI/Lx * (sphere(i,j)%lon - xc2)) yprime2 = 1.0D0/(DD_PI*sigmay_dbv) * SIN(DD_PI/Ly * (sphere(i,j)%lat - yc2)) p(i,j) = -g*dh_dbv * (EXP(-0.5D0 * (xprime1 * xprime1 + yprime1 * yprime1)) + EXP(-0.5D0 * (xprime2 * xprime2 + yprime2 * yprime2)) - 4.0D0*DD_PI*sigmax_dbv*sigmay_dbv) end do end do end function planar_dbl_vrtx_geopotential ! ====================================================== ! planar_dbl_vrtx_invariants: ! compute invariants on page 221 of Williamson, et. al. ! on a non-staggered grid. ! ====================================================== subroutine planar_dbl_vrtx_invariants(elem, iounit,tl,pmean,edge2,deriv,hybrid,nets,nete) type(element_t), intent(inout) :: elem(:) integer , intent(in) :: iounit type (TimeLevel_t) , intent(in) :: tl real (kind=real_kind), intent(in) :: pmean type (EdgeBuffer_t) :: edge2 type (derivative_t) :: deriv type (hybrid_t) , intent(in) :: hybrid integer , intent(in) :: nets,nete ! Local variables integer :: ie,i,j,kptr real (kind=real_kind) :: vco(np,np,2) real (kind=real_kind) :: gv(np,np,2) real (kind=real_kind) :: zeta(np,np,nets:nete) real (kind=real_kind) :: E(np,np) real (kind=real_kind) :: mass(np,np,nets:nete) real (kind=real_kind) :: energy(np,np,nets:nete) real (kind=real_kind) :: penst(np,np,nets:nete) real (kind=real_kind) :: vort(np,np,nets:nete) real (kind=real_kind) :: div(np,np,nets:nete) real (kind=real_kind) :: divsq(np,np,nets:nete) real (kind=real_kind) :: v1,v2 real (kind=real_kind) :: h,hstar,hs real (kind=real_kind) :: time real (kind=real_kind) :: Imass,Ienergy,Ipenst,Ivort,Idiv,Idivsq integer :: npts integer :: n0 n0 = tl%n0 if (tl%nstep == 0) then if (hybrid%masterthread) then open(iounit+0,file=TRIM(output_prefix)//"planar_dbl_vrtx.mass",status="unknown",form="formatted") open(iounit+1,file=TRIM(output_prefix)//"planar_dbl_vrtx.energy",status="unknown",form="formatted") open(iounit+2,file=TRIM(output_prefix)//"planar_dbl_vrtx.penst",status="unknown",form="formatted") open(iounit+3,file=TRIM(output_prefix)//"planar_dbl_vrtx.vort",status="unknown",form="formatted") open(iounit+4,file=TRIM(output_prefix)//"planar_dbl_vrtx.div",status="unknown",form="formatted") open(iounit+5,file=TRIM(output_prefix)//"planar_dbl_vrtx.rmsdiv",status="unknown",form="formatted") endif endif do ie=nets,nete do j=1,np do i=1,np v1 = elem(ie)%state%v(i,j,1,1,n0) v2 = elem(ie)%state%v(i,j,2,1,n0) vco(i,j,1) = elem(ie)%met(i,j,1,1)*v1 + elem(ie)%met(i,j,1,2)*v2 vco(i,j,2) = elem(ie)%met(i,j,2,1)*v1 + elem(ie)%met(i,j,2,2)*v2 gv(i,j,1) = elem(ie)%metdet(i,j)*v1 gv(i,j,2) = elem(ie)%metdet(i,j)*v2 end do end do zeta(:,:,ie) = vorticity(vco,deriv) div(:,:,ie) = divergence(gv,deriv) do j=1,np do i=1,np zeta(i,j,ie)=elem(ie)%mp(i,j)*(zeta(i,j,ie)/elem(ie)%metdet(i,j)) div(i,j,ie)=elem(ie)%mp(i,j)*(div(i,j,ie)/elem(ie)%metdet(i,j)) end do end do kptr=0 call edgeVpack(edge2, zeta(1,1,ie), 1, kptr,ie) kptr=1 call edgeVpack(edge2, div(1,1,ie), 1, kptr,ie) end do call bndry_exchangeV(hybrid,edge2) do ie=nets,nete kptr=0 call edgeVunpack(edge2, zeta(1,1,ie), 1, kptr, ie) kptr=1 call edgeVunpack(edge2, div(1,1,ie), 1, kptr, ie) do j=1,np do i=1,np zeta(i,j,ie) = elem(ie)%rmp(i,j)*zeta(i,j,ie) div(i,j,ie) = elem(ie)%rmp(i,j)*div(i,j,ie) penst(i,j,ie) = 0.5D0*(zeta(i,j,ie) + elem(ie)%fcor(i,j))**2 end do end do do j=1,np do i=1,np v1 = elem(ie)%state%v(i,j,1,1,n0) v2 = elem(ie)%state%v(i,j,2,1,n0) vco(i,j,1) = elem(ie)%met(i,j,1,1)*v1 + elem(ie)%met(i,j,1,2)*v2 vco(i,j,2) = elem(ie)%met(i,j,2,1)*v1 + elem(ie)%met(i,j,2,2)*v2 E(i,j) = 0.5D0*( vco(i,j,1)*v1 + vco(i,j,2)*v2 ) end do end do do j=1,np do i=1,np hstar = (elem(ie)%state%p(i,j,1,n0) + pmean)/g hs = elem(ie)%state%ps(i,j)/g h = hstar + hs mass(i,j,ie) = hstar energy(i,j,ie) = hstar*E(i,j) + 0.5D0*g*(h**2 - hs**2) penst(i,j,ie) = penst(i,j,ie)/hstar divsq(i,j,ie) = div(i,j,ie)*div(i,j,ie) vort(i,j,ie) = zeta(i,j,ie) end do end do end do npts=np Imass = global_integral(elem,mass(:,:,nets:nete),hybrid,npts,nets,nete) Ienergy = global_integral(elem,energy(:,:,nets:nete),hybrid,npts,nets,nete) Ipenst = global_integral(elem,penst(:,:,nets:nete),hybrid,npts,nets,nete) Ivort = global_integral(elem,vort(:,:,nets:nete),hybrid,npts,nets,nete) Idiv = global_integral(elem,div(:,:,nets:nete),hybrid,npts,nets,nete) Idivsq = global_integral(elem,divsq(:,:,nets:nete),hybrid,npts,nets,nete) time = Time_at(tl%nstep) if (hybrid%masterthread) then !print *,'mass, pmean/g = ',Imass,pmean/g if (tl%nstep==0) then Imass0 = Imass Ienergy0 = Ienergy Ipenst0 = Ipenst end if write(iounit+0,*)time/secpday,(Imass-Imass0)/Imass0,Imass0 write(iounit+1,*)time/secpday,(Ienergy-Ienergy0)/Ienergy0,Ienergy0 write(iounit+2,*)time/secpday,(Ipenst-Ipenst0)/Ipenst0,Ipenst0 write(iounit+3,*)time/secpday,Ivort write(iounit+4,*)time/secpday,Idiv write(iounit+5,*)time/secpday,SQRT(Idivsq) write(6,*)"" write(6,*)time/secpday,"mass =",(Imass-Imass0)/Imass0 write(6,*)time/secpday,"total energy =",(Ienergy-Ienergy0)/Ienergy0 write(6,*)time/secpday,"pot enstrophy =",(Ipenst-Ipenst0)/Ipenst0 write(6,*)time/secpday,"vorticity =",Ivort write(6,*)time/secpday,"divergence =",Idiv write(6,*)time/secpday,"RMS divergence=",SQRT(Idivsq) end if end subroutine planar_dbl_vrtx_invariants end module shallow_water_mod