diff --git a/doc/getm_pub.bib b/doc/getm_pub.bib index 1d27b0f2..47516467 100644 --- a/doc/getm_pub.bib +++ b/doc/getm_pub.bib @@ -42,6 +42,16 @@ doi = {10.1029/2020JC016411}, } +@Article{fofonova2021pst, + author = {Fofonova, Vera and K\"arn\"a, Tuomas and Klingbeil, Knut and Androsov, Alexey and Kuznetsov, Ivan and Sidorenko, Dmitry and Danilov, Sergey and Burchard, Hans and Wiltshire, Karen Helen}, + journal = {Geoscientific Model Development}, + title = {{Plume spreading test case for coastal ocean models}}, + year = {2021}, + pages = {6945--6975}, + volume = {14}, + doi = {10.5194/gmd-14-6945-2021}, +} + @article{Friedland2021, author = {Friedland, Ren{\'{e}} and Macias, Diego and Cossarini, Gianpiero and Daewel, Ute and Estournel, Claude and Garcia-Gorriz, Elisa and Gr{\'{e}}goire, Marilaure and Gustafson, Bo and Grizzetti, Bruna and Kalaroni, Sofia and Kerimoglu, Onur and Lazzari, Paolo and Lenhart, Hermann and Lessin, Gennadi and Maljutenko, Ilja and Miladinova, Svetla and M{\"{u}}ller-Karulis, B{\"{a}}rbel and Neumann, Thomas and Parn, Ove and P{\"{a}}tsch, Johannes and Piroddi, Chiara and Raudsepp, Urmas and Schrum, Corinna and Stegert, Christoph and Stips, Adolf and Tsiaras, Kostas and Ulses, Caroline and Vandenbulcke, Luc}, doi = {10.3389/fmars.2021.596126}, diff --git a/schemas/getm-2.5.defaults b/schemas/getm-2.5.defaults index 1346da0a..6311a3ae 100644 --- a/schemas/getm-2.5.defaults +++ b/schemas/getm-2.5.defaults @@ -620,6 +620,15 @@ -1.0 + + 0 + + + 0.0 + + + 40.0 + 0 diff --git a/schemas/getm-2.5.schema b/schemas/getm-2.5.schema index bacc0c8b..edf2df51 100644 --- a/schemas/getm-2.5.schema +++ b/schemas/getm-2.5.schema @@ -990,6 +990,9 @@ + + + diff --git a/src/2d/adv_split_u.F90 b/src/2d/adv_split_u.F90 index d9caf7d3..14bbb4ab 100644 --- a/src/2d/adv_split_u.F90 +++ b/src/2d/adv_split_u.F90 @@ -151,6 +151,7 @@ ! !USES: use domain, only: imin,imax,jmin,jmax use util, only: adv_reconstruct + use advection, only: ffluxu2 !$ use omp_lib IMPLICIT NONE ! @@ -273,6 +274,7 @@ if (associated(ffluxu)) then ffluxu(_U2DFIELD_) = ffluxu(_U2DFIELD_) + splitfac*uflux(_U2DFIELD_) end if + ffluxu2(_U2DFIELD_) = ffluxu2(_U2DFIELD_) + splitfac*uflux2(_U2DFIELD_) fio = fi Dio = Di !$OMP END WORKSHARE diff --git a/src/2d/adv_split_v.F90 b/src/2d/adv_split_v.F90 index cc73fa62..cc072789 100644 --- a/src/2d/adv_split_v.F90 +++ b/src/2d/adv_split_v.F90 @@ -21,6 +21,7 @@ ! !USES: use domain, only: imin,imax,jmin,jmax use util, only: adv_reconstruct + use advection, only: ffluxv2 !$ use omp_lib IMPLICIT NONE ! @@ -127,6 +128,7 @@ if (associated(ffluxv)) then ffluxv(_V2DFIELD_) = ffluxv(_V2DFIELD_) + splitfac*vflux(_V2DFIELD_) end if + ffluxv2(_V2DFIELD_) = ffluxv2(_V2DFIELD_) + splitfac*vflux2(_V2DFIELD_) fio = fi Dio = Di !$OMP END WORKSHARE diff --git a/src/2d/advection.F90 b/src/2d/advection.F90 index 024c1594..e97854a1 100644 --- a/src/2d/advection.F90 +++ b/src/2d/advection.F90 @@ -82,10 +82,12 @@ logical,dimension(E2DFIELD),target :: mask_updateH logical,dimension(E2DFIELD),target :: mask_uflux,mask_vflux,mask_xflux logical,dimension(E2DFIELD),target :: mask_uupdateU,mask_vupdateV + REALTYPE,dimension(E2DFIELD),public :: ffluxu2,ffluxv2 #else logical,dimension(:,:),allocatable,target :: mask_updateH logical,dimension(:,:),allocatable,target :: mask_uflux,mask_vflux,mask_xflux logical,dimension(:,:),allocatable,target :: mask_uupdateU,mask_vupdateV + REALTYPE,dimension(:,:),allocatable,public :: ffluxu2,ffluxv2 #endif ! ! !REVISION HISTORY: @@ -234,6 +236,12 @@ allocate(mask_vupdateV(E2DFIELD),stat=rc) ! work array if (rc /= 0) stop 'init_advection: Error allocating memory (mask_vupdateV)' + + allocate(ffluxu2(E2DFIELD),stat=rc) ! work array + if (rc /= 0) stop 'init_advection: Error allocating memory (ffluxu2)' + + allocate(ffluxv2(E2DFIELD),stat=rc) ! work array + if (rc /= 0) stop 'init_advection: Error allocating memory (ffluxv2)' #endif ! Note (KK): In this module pointers are used extensively. @@ -468,6 +476,8 @@ end if end if + ffluxu2 = _ZERO_ ; ffluxv2 = _ZERO_ + if (present(Dires)) then p_Di => Dires else diff --git a/src/2d/variables_2d.F90 b/src/2d/variables_2d.F90 index f63d4585..08b0cb60 100644 --- a/src/2d/variables_2d.F90 +++ b/src/2d/variables_2d.F90 @@ -34,6 +34,7 @@ logical :: calc_taubmax=.false. logical :: save_numdis_2d=.false. logical :: save_phydis_2d=.false. + logical :: save_rvol=.false. #ifdef STATIC #include "static_2d.h" @@ -61,6 +62,8 @@ REALTYPE,dimension(:,:),contiguous,pointer :: numdis_2d=>null() REALTYPE,dimension(:,:),contiguous,pointer :: phydis_2d=>null() + REALTYPE,dimension(:,:),allocatable :: rvol + integer :: size2d_field integer :: mem2d ! @@ -317,6 +320,12 @@ taubmax = _ZERO_ end if + if (save_rvol) then + allocate(rvol(E2DFIELD),stat=rc) + if (rc /= 0) stop 'postinit_2d: Error allocating memory (rvol)' + rvol = _ZERO_ + end if + end if #ifdef DEBUG @@ -340,6 +349,7 @@ ! !USES: use field_manager use meteo, only: fwf_method + use domain, only: nriverl IMPLICIT NONE ! ! !INPUT PARAMETERS: @@ -412,6 +422,10 @@ call fm%register('fwf', 'm', 'surface freshwater fluxes', standard_name='', data2d=fwf(_2D_FO_), category='2d', add_to_mean_before_save=.true., output_level=output_level_debug) end if + if (nriverl .gt. 0) then + call fm%register('rvol', 'm**3', 'river discharge', category='2d', fill_value=_ZERO_, output_level=output_level_debug, used=save_rvol) + end if + return end subroutine register_2d_variables !EOC @@ -445,6 +459,7 @@ if (associated(taubmax)) call fm%send_data('taubmax', taubmax(_2D_FO_)) if (associated(numdis_2d)) call fm%send_data('numdis_2d', numdis_2d(_2D_FO_)) if (associated(phydis_2d)) call fm%send_data('phydis_2d', phydis_2d(_2D_FO_)) + if (allocated(rvol)) call fm%send_data('rvol', rvol(_2D_FO_)) return end subroutine finalize_register_2d_variables diff --git a/src/3d/advection_3d.F90 b/src/3d/advection_3d.F90 index 8e0f855c..f7d2f2e5 100644 --- a/src/3d/advection_3d.F90 +++ b/src/3d/advection_3d.F90 @@ -45,6 +45,12 @@ "full step splitting: u + v + w ", & "half step splitting: u/2 + v/2 + w + v/2 + u/2", & "hor/ver splitting: uv + w "/) + +#ifdef STATIC + REALTYPE,dimension(I3DFIELD),public :: ffluxu2_3d, ffluxv2_3d +#else + REALTYPE,dimension(:,:,:),allocatable,public :: ffluxu2_3d, ffluxv2_3d +#endif ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Hans Burchard @@ -86,6 +92,9 @@ ! !USES: IMPLICIT NONE ! +! !LOCAL VARIABLES: + integer :: rc +! !EOP !------------------------------------------------------------------------- !BOC @@ -106,6 +115,14 @@ end if #endif +#ifndef STATIC + allocate(ffluxu2_3d(I3DFIELD),stat=rc) ! work array + if (rc /= 0) stop 'init_advection_3d: Error allocating memory (ffluxu2_3d)' + + allocate(ffluxv2_3d(I3DFIELD),stat=rc) ! work array + if (rc /= 0) stop 'init_advection_3d: Error allocating memory (ffluxv2_3d)' +#endif + #ifdef DEBUG write(debug,*) 'Leaving init_advection_3d()' write(debug,*) @@ -376,6 +393,8 @@ end if end if + ffluxu2_3d = _ZERO_ ; ffluxv2_3d = _ZERO_ + if (present(hires)) then p_hi => hires else @@ -405,12 +424,15 @@ if (calc_ffluxu) pa_ffluxu2d(k)%p2d = _ZERO_ if (calc_ffluxv) pa_ffluxv2d(k)%p2d = _ZERO_ #endif + ffluxu2 = _ZERO_ ; ffluxv2 = _ZERO_ call do_advection(dt,f(:,:,k),uu(:,:,k),vv(:,:,k), & hu(:,:,k),hv(:,:,k),ho(:,:,k),hn(:,:,k), & split,hscheme,AH,tag2d, & Dires=p_hi(:,:,k),advres=p_adv(:,:,k), & ffluxu=pa_ffluxu2d(k)%p2d,ffluxv=pa_ffluxv2d(k)%p2d, & nvd=pa_nvd2d(k)%p2d) + ffluxu2_3d(:,:,k) = ffluxu2 + ffluxv2_3d(:,:,k) = ffluxv2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxu) ffluxu(:,:,k) = pa_ffluxu2d(k)%p2d @@ -435,19 +457,23 @@ if (calc_ffluxu) pa_ffluxu2d(k)%p2d = _ZERO_ if (calc_ffluxv) pa_ffluxv2d(k)%p2d = _ZERO_ #endif + ffluxu2 = _ZERO_ call adv_split_u(dt,f(:,:,k),fi(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),uu(:,:,k),hu(:,:,k), & adv_grid%dxu,adv_grid%dyu,adv_grid%arcd1, & _ONE_,hscheme,AH, & adv_grid%mask_uflux,adv_grid%mask_uupdate, & pa_ffluxu2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxu2_3d(:,:,k) = ffluxu2 #ifndef SLICE_MODEL + ffluxv2 = _ZERO_ call adv_split_v(dt,f(:,:,k),fi(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),vv(:,:,k),hv(:,:,k), & adv_grid%dxv,adv_grid%dyv,adv_grid%arcd1, & _ONE_,hscheme,AH, & adv_grid%mask_vflux,adv_grid%mask_vupdate, & pa_ffluxv2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxv2_3d(:,:,k) = ffluxv2 #endif #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d @@ -532,12 +558,15 @@ if (calc_ffluxu) pa_ffluxu2d(k)%p2d = _ZERO_ if (calc_ffluxv) pa_ffluxv2d(k)%p2d = _ZERO_ #endif + ffluxu2 = _ZERO_ ; ffluxv2 = _ZERO_ call do_advection(dt,f(:,:,k),uu(:,:,k),vv(:,:,k), & hu(:,:,k),hv(:,:,k),ho(:,:,k),hn(:,:,k), & FULLSPLIT,hscheme,AH,tag2d, & Dires=p_hi(:,:,k),advres=p_adv(:,:,k), & ffluxu=pa_ffluxu2d(k)%p2d,ffluxv=pa_ffluxv2d(k)%p2d, & nvd=pa_nvd2d(k)%p2d) + ffluxu2_3d(:,:,k) = ffluxu2 + ffluxv2_3d(:,:,k) = ffluxv2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxu) ffluxu(:,:,k) = pa_ffluxu2d(k)%p2d @@ -560,12 +589,14 @@ if (calc_nvd) pa_nvd2d (k)%p2d = _ZERO_ if (calc_ffluxu) pa_ffluxu2d(k)%p2d = _ZERO_ #endif + ffluxu2 = _ZERO_ call adv_split_u(dt,f(:,:,k),f(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),uu(:,:,k),hu(:,:,k), & adv_grid%dxu,adv_grid%dyu,adv_grid%arcd1, & _HALF_,hscheme,AH, & adv_grid%mask_uflux,adv_grid%mask_uupdate, & pa_ffluxu2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxu2_3d(:,:,k) = ffluxu2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxu) ffluxu(:,:,k) = pa_ffluxu2d(k)%p2d @@ -585,12 +616,14 @@ if (calc_nvd) pa_nvd2d (k)%p2d = nvd(:,:,k) if (calc_ffluxv) pa_ffluxv2d(k)%p2d = _ZERO_ #endif + ffluxv2 = _ZERO_ call adv_split_v(dt,f(:,:,k),f(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),vv(:,:,k),hv(:,:,k), & adv_grid%dxv,adv_grid%dyv,adv_grid%arcd1, & _HALF_,hscheme,AH, & adv_grid%mask_vflux,adv_grid%mask_vupdate, & pa_ffluxv2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxv2_3d(:,:,k) = ffluxv2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxv) ffluxv(:,:,k) = pa_ffluxv2d(k)%p2d @@ -616,12 +649,14 @@ if (calc_nvd) pa_nvd2d (k)%p2d = nvd (:,:,k) if (calc_ffluxv) pa_ffluxv2d(k)%p2d = ffluxv(:,:,k) #endif + ffluxv2 = ffluxv2_3d(:,:,k) call adv_split_v(dt,f(:,:,k),f(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),vv(:,:,k),hv(:,:,k), & adv_grid%dxv,adv_grid%dyv,adv_grid%arcd1, & _HALF_,hscheme,AH, & adv_grid%mask_vflux,adv_grid%mask_vupdate, & pa_ffluxv2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxv2_3d(:,:,k) = ffluxv2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxv) ffluxv(:,:,k) = pa_ffluxv2d(k)%p2d @@ -641,12 +676,14 @@ if (calc_nvd) pa_nvd2d (k)%p2d = nvd (:,:,k) if (calc_ffluxu) pa_ffluxu2d(k)%p2d = ffluxu(:,:,k) #endif + ffluxu2 = ffluxu2_3d(:,:,k) call adv_split_u(dt,f(:,:,k),f(:,:,k),p_hi(:,:,k), & p_adv(:,:,k),uu(:,:,k),hu(:,:,k), & adv_grid%dxu,adv_grid%dyu,adv_grid%arcd1, & _HALF_,hscheme,AH, & adv_grid%mask_uflux,adv_grid%mask_uupdate, & pa_ffluxu2d(k)%p2d,pa_nvd2d(k)%p2d) + ffluxu2_3d(:,:,k) = ffluxu2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxu) ffluxu(:,:,k) = pa_ffluxu2d(k)%p2d @@ -671,13 +708,15 @@ if (calc_ffluxu) pa_ffluxu2d(k)%p2d = _ZERO_ if (calc_ffluxv) pa_ffluxv2d(k)%p2d = _ZERO_ #endif + ffluxu2 = _ZERO_ ; ffluxv2 = _ZERO_ call do_advection(dt,f(:,:,k),uu(:,:,k),vv(:,:,k), & hu(:,:,k),hv(:,:,k),ho(:,:,k),hn(:,:,k), & NOSPLIT,hscheme,AH,tag2d, & Dires=p_hi(:,:,k),advres=p_adv(:,:,k), & ffluxu=pa_ffluxu2d(k)%p2d,ffluxv=pa_ffluxv2d(k)%p2d, & nvd=pa_nvd2d(k)%p2d) - + ffluxu2_3d(:,:,k) = ffluxu2 + ffluxv2_3d(:,:,k) = ffluxv2 #ifndef _POINTER_REMAP_ if (calc_nvd) nvd (:,:,k) = pa_nvd2d (k)%p2d if (calc_ffluxu) ffluxu(:,:,k) = pa_ffluxu2d(k)%p2d diff --git a/src/3d/getm_fabm.F90 b/src/3d/getm_fabm.F90 index a3859072..408760c6 100644 --- a/src/3d/getm_fabm.F90 +++ b/src/3d/getm_fabm.F90 @@ -70,7 +70,7 @@ type t_pa3d REALTYPE,dimension(:,:,:),pointer,contiguous :: p3d end type t_pa3d - type(t_pa3d),dimension(:),allocatable :: pa_nummix_fabm_pel,pa_phymix_fabm_pel + type(t_pa3d),dimension(:),allocatable,public :: pa_nummix_fabm_pel,pa_phymix_fabm_pel integer :: fabm_adv_split=0 integer :: fabm_adv_hor=1 integer :: fabm_adv_ver=1 diff --git a/src/3d/m3d.F90 b/src/3d/m3d.F90 index e80361cb..83da4a83 100644 --- a/src/3d/m3d.F90 +++ b/src/3d/m3d.F90 @@ -329,7 +329,7 @@ else T = _ZERO_ ; S = _ZERO_ ; rho = _ZERO_ if(calc_temp) call init_temperature(hotstart) - if(calc_salt) call init_salinity(hotstart) + if(calc_salt) call init_salinity(runtype,hotstart) if (runtype .eq. 4) then LEVEL2 'AH_min_method = ', AH_min_method select case (AH_min_method) diff --git a/src/3d/rivers.F90 b/src/3d/rivers.F90 index f12e27b7..080a610c 100644 --- a/src/3d/rivers.F90 +++ b/src/3d/rivers.F90 @@ -28,11 +28,13 @@ ! !USES: use domain, only: imin,jmin,imax,jmax,ioff,joff use domain, only: H,az,kmax,arcd1 + use domain, only: nriverl,ri,rj use time, only: write_time_string,timestr use variables_2d, only: dtm,z #ifndef NO_BAROCLINIC - use m3d, only: update_salt,update_temp + use m3d, only: update_salt,update_temp,dt use variables_3d, only: hn,ssen,T,S + use variables_3d, only: phymix_S #endif #ifdef GETM_BIO use bio, only: bio_calc @@ -76,6 +78,7 @@ #ifdef _FABM_ REALTYPE, public, allocatable :: river_fabm(:,:) #endif + integer, public, allocatable :: rnl(:) ! ! !PRIVATE DATA MEMBERS: integer :: river_format=2 @@ -222,9 +225,11 @@ else if (.not. outside) then xxx = ' inside' ok(n) = 1 + nriverl = nriverl + 1 else xxx = ' in halo' ok(n) = 2 + nriverl = nriverl + 1 end if bathy = H(i,j) if (rzu(n) .gt. rzl(n)) then @@ -249,6 +254,23 @@ LEVEL3 trim(line) end do + allocate(ri(nriverl),stat=rc) + if (rc /= 0) stop 'rivers: Error allocating memory (ri)' + allocate(rj(nriverl),stat=rc) + if (rc /= 0) stop 'rivers: Error allocating memory (rj)' + allocate(rnl(nriverl),stat=rc) + if (rc /= 0) stop 'rivers: Error allocating memory (rnl)' + + nn = 0 + do n=1,nriver + if (ok(n).eq.1 .or. ok(n).eq.2) then + nn = nn + 1 + ri(nn) = ir(n) - ioff + rj(nn) = jr(n) - joff + rnl(nn) = n + end if + end do + ! Calculate the number of used gridboxes. ! This particular section is prepared for multi-cell rivers not-in-sequence. LEVEL3 'Number of unique rivers: ',rriver @@ -511,6 +533,9 @@ ! the layer heights are increased proportionally. ! ! !USES: + use variables_2d, only: rvol_2d=>rvol + use variables_3d, only: saltbins,Si_s,rvol_s + use getm_timers , only: tic,toc,TIM_TEF IMPLICIT NONE ! ! !INPUT PARAMETERS: @@ -518,10 +543,12 @@ logical, intent(in) :: do_3d ! ! !LOCAL VARIABLES: - integer :: i,j,k,m,n + integer :: i,j,k,l,m,n,nn integer :: kl,kh REALTYPE :: rvol,height REALTYPE :: river_depth,x + REALTYPE :: S_old(0:kmax),S_riv(0:kmax) + integer :: binidx(0:kmax) !EOP !----------------------------------------------------------------------- !BOC @@ -548,16 +575,18 @@ select case (river_method) case(0) case(1,2) + nn = 0 do n=1,nriver if(ok(n) .gt. 0) then + nn = nn + 1 i = ir(n)-ioff; j = jr(n)-joff rvol = ramp * dtm * river_flow(n) * flow_fraction(n) + if (allocated(rvol_2d)) rvol_2d(i,j) = rvol if (ok(n) .eq. 1) then irr(n) = irr(n) + rvol end if height = rvol * ARCD1 z(i,j) = z(i,j) + height -#ifndef NO_BAROCLINIC macro_height(n)=macro_height(n)+height ! on macrotime step adjust 3d fields if (do_3d) then @@ -582,16 +611,40 @@ kh = m end if river_depth = sum(hn(i,j,kl:kh)) + +! Changes of total and layer height due to river inflow: + ssen(i,j) = ssen(i,j)+macro_height(n) + x = macro_height(n) / river_depth + hn(i,j,kl:kh) = hn(i,j,kl:kh) * ( _ONE_ + x ) + +#ifndef NO_BAROCLINIC if (macro_height(n).gt._ZERO_ .or. .not.river_outflow_properties_follow_source_cell) then if (update_salt ) then + S_old(:) = S(i,j,:) if ( river_salt(n) .ne. salt_missing) then - S(i,j,kl:kh) = (S(i,j,kl:kh)*river_depth & - + river_salt(n)*macro_height(n)) & - / (river_depth+macro_height(n)) + S(i,j,kl:kh) = ( S(i,j,kl:kh) + x*river_salt(n) ) / ( _ONE_ + x ) + S_riv(kl:kh) = river_salt(n) else - S(i,j,kl:kh) = S(i,j,kl:kh)*river_depth & - / (river_depth+macro_height(n)) + S(i,j,kl:kh) = S(i,j,kl:kh) / ( _ONE_ + x ) + S_riv(kl:kh) = _ZERO_ end if +! only need to calculate mixing here (no mixing for zero-gradient) + if (associated(phymix_S)) then + phymix_S(i,j,kl:kh) = & + hn(i,j,kl:kh)*( ( S_old(kl:kh)**2 + x*S_riv(kl:kh)**2 ) / ( _ONE_ + x ) - S(i,j,kl:kh)**2 ) / dt + end if + if (allocated(rvol_s)) then + call tic(TIM_TEF) + rvol_s(i,j,:) = _ZERO_ + call rebin(kmax,S_old,saltbins,Si_s,binidx) + do k=kl,kh + l = binidx(k) + if (l .lt. 0) cycle + rvol_s(i,j,l) = rvol_s(i,j,l) + x*hn(i,j,k)/((_ONE_+x)*ARCD1) + end do + call toc(TIM_TEF) + end if + end if if (update_temp .and. river_temp(n) .ne. temp_missing) then T(i,j,kl:kh) = (T(i,j,kl:kh)*river_depth & @@ -623,14 +676,11 @@ end if #endif end if +#endif -! Changes of total and layer height due to river inflow: - hn(i,j,kl:kh) = hn(i,j,kl:kh)/river_depth & - *(river_depth+macro_height(n)) - ssen(i,j) = ssen(i,j)+macro_height(n) macro_height(n) = _ZERO_ end if -#endif + end if end do case default diff --git a/src/3d/salinity.F90 b/src/3d/salinity.F90 index 99c7a9ae..eadd437e 100644 --- a/src/3d/salinity.F90 +++ b/src/3d/salinity.F90 @@ -20,6 +20,7 @@ !KB use get_field, only: get_3d_field use variables_2d, only: fwf_int use variables_3d, only: rk,S,hn,kmin,sss + use variables_3d, only: saltbins,saltbinmin,saltbinmax,Si_s use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG IMPLICIT NONE @@ -61,6 +62,7 @@ integer :: salt_check=0 REALTYPE :: min_salt=_ZERO_,max_salt=40*_ONE_ integer,public :: nonnegsalt_method=0 + integer, allocatable :: binidx_kij(:,:,:) ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Hans Burchard @@ -77,7 +79,7 @@ ! \label{sec-init-salinity} ! ! !INTERFACE: - subroutine init_salinity(hotstart) + subroutine init_salinity(runtype,hotstart) ! ! !DESCRIPTION: ! @@ -104,10 +106,11 @@ IMPLICIT NONE ! ! !INPUT PARAMETERS: + integer,intent(in) :: runtype logical,intent(in) :: hotstart ! ! !LOCAL VARIABLES: - integer :: i,j,k,n,rc + integer :: i,j,k,l,n,rc integer :: status NAMELIST /salt/ & salt_method,salt_const,salt_file, & @@ -115,6 +118,7 @@ salt_adv_split,salt_adv_hor,salt_adv_ver, & avmols,salt_AH_method,salt_AH_const,salt_AH_Prt, & salt_AH_stirr_const,sss_nudging_time, & + saltbins,saltbinmin,saltbinmax, & salt_check,min_salt,max_salt, & nonnegsalt_method !EOP @@ -136,10 +140,6 @@ LEVEL3 'avmols = ',real(avmols) end if - if (.not. hotstart) then - call init_salinity_field() - end if - ! Sanity checks for advection specifications LEVEL3 'Advection of salinity' if (salt_adv_hor .eq. J7) stop 'init_salinity: J7 not implemented yet' @@ -188,6 +188,22 @@ allocate(sss(E2DFIELD),stat=rc) end if + if (runtype.eq.4 .and. saltbins.gt.0) then + LEVEL3 'salinity space for TEF:',real(saltbinmin),real(saltbinmax),'(',saltbins,'bins)' + allocate(Si_s(0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (Si_s)' + do l=0,saltbins + Si_s(l) = saltbinmin + (saltbinmax-saltbinmin)/saltbins * l + end do + allocate(binidx_kij(_KRANGE_,I2DFIELD),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (binidx_kij)' + binidx_kij = -1 + end if + + if (.not. hotstart) then + call init_salinity_field() + end if + LEVEL3 'salt_check=',salt_check if (salt_check .ne. 0) then LEVEL4 'doing sanity check on salinity' @@ -298,6 +314,14 @@ call wait_halo(D_TAG) call mirror_bdy_3d(S,D_TAG) + if (allocated(binidx_kij)) then + do j = jmin-HALO,jmax+HALO + do i = imin-HALO,imax+HALO + if (az(i,j).ne.0) call rebin(kmax,S(i,j,:),saltbins,Si_s,binidx_kij(:,i,j)) + end do + end do + end if + #ifdef DEBUG write(debug,*) 'Leaving init_salinity_field()' write(debug,*) @@ -312,7 +336,7 @@ ! !IROUTINE: do_salinity - salinity equation \label{sec-do-salinity} ! ! !INTERFACE: - subroutine do_salinity(n) + subroutine do_salinity(loop) ! ! !DESCRIPTION: ! @@ -333,24 +357,40 @@ ! ! !USES: use advection_3d, only: do_advection_3d + use advection_3d, only: ffluxu2_3d,ffluxv2_3d use variables_3d, only: dt,cnpar,hn,ho,nuh,uu,vv,ww,hun,hvn - use domain, only: imin,imax,jmin,jmax,kmax,az - use getm_timers, only: tic,toc,TIM_SALT,TIM_SALTH,TIM_MIXANALYSIS + use domain, only: imin,imax,jmin,jmax,kmax,az,au,av + use domain, only: dyu,dxv + use domain, only: nriverl,ri,rj + use getm_timers, only: tic,toc,TIM_SALT,TIM_SALTH,TIM_MIXANALYSIS,TIM_TEF use variables_3d, only: Sfluxu,Sfluxv,Sfluxw - use variables_3d, only: nummix_S, phymix_S + use variables_3d, only: Sfluxu2,Sfluxv2 + use variables_3d, only: counts_s,flags_s,h_s,hS_s,hS2_s,hpmS_s,hnmS_s + use variables_3d, only: hu_s,uu_s,Sfluxu_s,S2fluxu_s + use variables_3d, only: hv_s,vv_s,Sfluxv_s,S2fluxv_s + use variables_3d, only: fwf_s,fwfS2_s, phymix_S_fwf_s, phymix_S_riv_s + use variables_3d, only: nummix_S, phymix_S, phymix_S_fwf !$ use omp_lib IMPLICIT NONE ! ! !INPUT PARAMETERS: - integer, intent(in) :: n + integer, intent(in) :: loop ! ! !LOCAL VARIABLES: - integer :: i,j,k,rc + integer :: i,j,k,l,n,rc REALTYPE, POINTER :: Res(:) REALTYPE, POINTER :: auxn(:),auxo(:) REALTYPE, POINTER :: a1(:),a2(:),a3(:),a4(:) REALTYPE, pointer :: fluxw(:) + REALTYPE :: sss_old(I2DFIELD) + REALTYPE, dimension(I3DFIELD) :: Sold,Sfluxu_tmp,Sfluxv_tmp + REALTYPE, dimension(I2DFIELD) :: wrk2d + REALTYPE :: rivmix(0:kmax,1:nriverl) + REALTYPE :: deltaS,hu_k,uu_k + REALTYPE :: S_k(0:kmax) + integer :: binidx(0:kmax) integer :: status + logical :: calc_phymix_fwf !EOP !----------------------------------------------------------------------- !BOC @@ -361,6 +401,15 @@ #endif call tic(TIM_SALT) + if (allocated(phymix_S_riv_s)) then + do n=1,nriverl + rivmix(:,n) = phymix_S(ri(n),rj(n),:) + end do + end if + + calc_phymix_fwf = ( associated(phymix_S) .or. allocated(phymix_S_fwf) ) + if (calc_phymix_fwf) sss_old = S(:,:,kmax) + do j=jmin-HALO,jmax+HALO do i=imin-HALO,imax+HALO if (az(i,j) .eq. 1) then @@ -369,28 +418,202 @@ end do end do + if (calc_phymix_fwf) then +! Note (KK): division by hn because we already refer to new volume +! where diffusion will be added + where (az .eq. 1) wrk2d = sss_old * S(:,:,kmax) * fwf_int / dt + if (associated(phymix_S)) then + where (az .eq. 1) phymix_S(:,:,kmax) = phymix_S(:,:,kmax) + wrk2d + end if + if (allocated(phymix_S_fwf)) then + where (az .eq. 1) phymix_S_fwf = wrk2d + end if + end if + + if (allocated(fwf_s).or.allocated(fwfS2_s)) then + call tic(TIM_TEF) + if (allocated(fwf_s )) fwf_s = _ZERO_ + if (allocated(fwfS2_s)) fwfS2_s = _ZERO_ + do j = jmin,jmax + do i = imin,imax + if (az(i,j) .eq. 1) then + l = binidx_kij(kmax,i,j) ! binning based on old salinity + if (l .ge. 0) then + if (allocated(fwf_s )) fwf_s (i,j,l) = fwf_int(i,j) + if (allocated(fwfS2_s)) fwfS2_s(i,j,l) = fwf_int(i,j)*S(i,j,kmax)*S(i,j,kmax) + end if + end if + end do + end do + call toc(TIM_TEF) + end if + + call do_advection_3d(dt,S,uu,vv,ww,hun,hvn,ho,hn, & salt_adv_split,salt_adv_hor,salt_adv_ver,_ZERO_,H_TAG, & ffluxu=Sfluxu,ffluxv=Sfluxv,ffluxw=Sfluxw,nvd=nummix_S) + if (allocated(Sfluxu2)) Sfluxu2 = ffluxu2_3d + if (allocated(Sfluxv2)) Sfluxv2 = ffluxv2_3d + +! KK-TODO: do we have to treat salt fluxes from different +! directional-split steps separately? + if (allocated(hu_s).or.allocated(uu_s).or.allocated(Sfluxu_s).or.allocated(S2fluxu_s)) then + call tic(TIM_TEF) + if (allocated(hu_s )) hu_s = _ZERO_ + if (allocated(uu_s )) uu_s = _ZERO_ + if (allocated(Sfluxu_s )) Sfluxu_s = _ZERO_ + if (allocated(S2fluxu_s)) S2fluxu_s = _ZERO_ + do j = jmin,jmax + do i = imin,imax + if (au(i,j).eq.1 .or. au(i,j).eq.2) then + where (uu(i,j,1:kmax) .ne. _ZERO_) + S_k(1:kmax) = Sfluxu(i,j,1:kmax) / uu(i,j,1:kmax) / DYU + else where + S_k(1:kmax) = saltbinmax + 9999.0_rk + end where + call rebin(kmax,S_k,saltbins,Si_s,binidx) + do k=1,kmax + l = binidx(k) + if (l .lt. 0) cycle + if (allocated(hu_s )) hu_s (i,j,l) = hu_s (i,j,l) + hun(i,j,k) + if (allocated(uu_s )) uu_s (i,j,l) = uu_s (i,j,l) + uu(i,j,k) + if (allocated(Sfluxu_s )) Sfluxu_s (i,j,l) = Sfluxu_s (i,j,l) + Sfluxu(i,j,k) + if (allocated(S2fluxu_s)) S2fluxu_s(i,j,l) = S2fluxu_s(i,j,l) + Sfluxu2(i,j,k) + end do + end if + end do + end do + call toc(TIM_TEF) + end if + if (allocated(hv_s).or.allocated(vv_s).or.allocated(Sfluxv_s).or.allocated(S2fluxv_s)) then + call tic(TIM_TEF) + if (allocated(hv_s )) hv_s = _ZERO_ + if (allocated(vv_s )) vv_s = _ZERO_ + if (allocated(Sfluxv_s )) Sfluxv_s = _ZERO_ + if (allocated(S2fluxv_s)) S2fluxv_s = _ZERO_ + do j = jmin,jmax + do i = imin,imax + if (av(i,j).eq.1 .or. av(i,j).eq.2) then + where (vv(i,j,1:kmax) .ne. _ZERO_) + S_k(1:kmax) = Sfluxv(i,j,1:kmax) / vv(i,j,1:kmax) / DXV + else where + S_k(1:kmax) = saltbinmax + 9999.0_rk + end where + call rebin(kmax,S_k,saltbins,Si_s,binidx) + do k=1,kmax + l = binidx(k) + if (l .lt. 0) cycle + if (allocated(hv_s )) hv_s (i,j,l) = hv_s (i,j,l) + hvn(i,j,k) + if (allocated(vv_s )) vv_s (i,j,l) = vv_s (i,j,l) + vv(i,j,k) + if (allocated(Sfluxv_s )) Sfluxv_s (i,j,l) = Sfluxv_s (i,j,l) + Sfluxv(i,j,k) + if (allocated(S2fluxv_s)) S2fluxv_s(i,j,l) = S2fluxv_s(i,j,l) + Sfluxv2(i,j,k) + end do + end if + end do + end do + call toc(TIM_TEF) + end if + if (salt_AH_method .gt. 0) then + ! S is not halo updated after advection call tic(TIM_SALTH) call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG) call wait_halo(D_TAG) call toc(TIM_SALTH) + if ( allocated(hu_s).or.allocated(uu_s).or.allocated(Sfluxu_s).or.allocated(S2fluxu_s) & + .or.allocated(hv_s).or.allocated(vv_s).or.allocated(Sfluxv_s).or.allocated(S2fluxv_s)) then + call tic(TIM_TEF) + Sold = S + do j = jmin,jmax+1 + do i = imin,imax+1 + if (az(i,j).ne.0) call rebin(kmax,S(i,j,:),saltbins,Si_s,binidx_kij(:,i,j)) + end do + end do + call toc(TIM_TEF) + end if + if (allocated(hu_s).or.allocated(uu_s).or.allocated(Sfluxu_s)) Sfluxu_tmp = Sfluxu + if (allocated(hv_s).or.allocated(vv_s).or.allocated(Sfluxv_s)) Sfluxv_tmp = Sfluxv + call tracer_diffusion(S,hn,salt_AH_method,salt_AH_const,salt_AH_Prt,salt_AH_stirr_const, & ffluxu=Sfluxu,ffluxv=Sfluxv, & phymix=phymix_S) - end if + if (allocated(Sfluxu2)) Sfluxu2 = Sfluxu2 + ffluxu2_3d + if (allocated(Sfluxv2)) Sfluxv2 = Sfluxv2 + ffluxv2_3d + + if (allocated(hu_s).or.allocated(uu_s).or.allocated(Sfluxu_s).or.allocated(S2fluxu_s)) then + call tic(TIM_TEF) + Sfluxu_tmp = Sfluxu - Sfluxu_tmp + do j = jmin,jmax + do i = imin,imax + if (au(i,j).eq.1 .or. au(i,j).eq.2) then + do k=1,kmax + if (Sfluxu_tmp(i,j,k) .ne. _ZERO_) then + deltaS = Sold(i+1,j,k) - Sold(i,j,k) + !if (deltaS .ne. _ZERO_) then + hu_k = _HALF_ * ( hn(i,j,k) + hn(i+1,j,k) ) + uu_k = Sfluxu_tmp(i,j,k) / deltaS / DYU +! fluxes with S(i) + l = binidx_kij(k,i,j) + if (l .ge. 0) then + if (allocated(hu_s )) hu_s (i,j,l) = hu_s (i,j,l) + hu_k + if (allocated(uu_s )) uu_s (i,j,l) = uu_s (i,j,l) - uu_k + if (allocated(Sfluxu_s )) Sfluxu_s (i,j,l) = Sfluxu_s (i,j,l) - uu_k*DYU*Sold(i,j,k) + if (allocated(S2fluxu_s)) S2fluxu_s(i,j,l) = S2fluxu_s(i,j,l) - uu_k*DYU*Sold(i,j,k)*Sold(i,j,k) + end if +! fluxes with S(i+1) + l = binidx_kij(k,i+1,j) + if (l .ge. 0) then + if (allocated(hu_s )) hu_s (i,j,l) = hu_s (i,j,l) + hu_k + if (allocated(uu_s )) uu_s (i,j,l) = uu_s (i,j,l) + uu_k + if (allocated(Sfluxu_s )) Sfluxu_s (i,j,l) = Sfluxu_s (i,j,l) + uu_k*DYU*Sold(i+1,j,k) + if (allocated(S2fluxu_s)) S2fluxu_s(i,j,l) = S2fluxu_s(i,j,l) + uu_k*DYU*Sold(i+1,j,k)*Sold(i+1,j,k) + end if + end if + end do + end if + end do + end do + call toc(TIM_TEF) + end if + if (allocated(hv_s).or.allocated(vv_s).or.allocated(Sfluxv_s).or.allocated(S2fluxv_s)) then + call tic(TIM_TEF) + Sfluxv_tmp = Sfluxv - Sfluxv_tmp + do j = jmin,jmax + do i = imin,imax + if (av(i,j).eq.1 .or. av(i,j).eq.2) then + do k=1,kmax + if (Sfluxv_tmp(i,j,k) .ne. _ZERO_) then + deltaS = Sold(i,j+1,k) - Sold(i,j,k) + !if (deltaS .ne. _ZERO_) then + hu_k = _HALF_ * ( hn(i,j,k) + hn(i,j+1,k) ) + uu_k = Sfluxv_tmp(i,j,k) / deltaS / DXV +! fluxes with S(j) + l = binidx_kij(k,i,j) + if (l .ge. 0) then + if (allocated(hv_s )) hv_s (i,j,l) = hv_s (i,j,l) + hu_k + if (allocated(vv_s )) vv_s (i,j,l) = vv_s (i,j,l) - uu_k + if (allocated(Sfluxv_s )) Sfluxv_s (i,j,l) = Sfluxv_s (i,j,l) - uu_k*DXV*Sold(i,j,k) + if (allocated(S2fluxv_s)) S2fluxv_s(i,j,l) = S2fluxv_s(i,j,l) - uu_k*DXV*Sold(i,j,k)*Sold(i,j,k) + end if +! fluxes with S(j+1) + l = binidx_kij(k,i,j+1) + if (l .ge. 0) then + if (allocated(hv_s )) hv_s (i,j,l) = hv_s (i,j,l) + hu_k + if (allocated(vv_s )) vv_s (i,j,l) = vv_s (i,j,l) + uu_k + if (allocated(Sfluxv_s )) Sfluxv_s (i,j,l) = Sfluxv_s (i,j,l) + uu_k*DXV*Sold(i,j+1,k) + if (allocated(S2fluxv_s)) S2fluxv_s(i,j,l) = S2fluxv_s(i,j,l) + uu_k*DXV*Sold(i,j+1,k)*Sold(i,j+1,k) + end if + end if + end do + end if + end do + end do + call toc(TIM_TEF) + end if - if (associated(phymix_S)) then - call toc(TIM_SALT) - call tic(TIM_MIXANALYSIS) - call physical_mixing(S,avmols,phymix_S,salt_AH_method) - call toc(TIM_MIXANALYSIS) - call tic(TIM_SALT) end if @@ -417,7 +640,7 @@ if (rc /= 0) stop 'do_salinity: Error allocating memory (auxo)' fluxw => null() - if (associated(Sfluxw)) then + if (associated(Sfluxw) .or. associated(phymix_S)) then allocate(fluxw(0:kmax),stat=rc) ! work array if (rc /= 0) stop 'do_salinity: Error allocating memory (fluxw)' end if @@ -486,6 +709,28 @@ if (associated(Sfluxw)) then Sfluxw(i,j,1:kmax-1) = Sfluxw(i,j,1:kmax-1) + fluxw(1:kmax-1) end if + if (associated(phymix_S)) then +! KK-TODO: only consider changes of S2 due to diffusion! +! (no sss_nudging or radiation for temperature) + do k=1,kmax-1 +#if 0 + fluxw(k) = fluxw(k) * & + ( cnpar *(Res (k)+Res (k+1)) & + + (_ONE_-cnpar)*(S(i,j,k)+S(i,j,k+1))) +#else + fluxw(k) = - ( & + auxo(k) * ( S(i,j,k+1) - S(i,j,k) ) * (S(i,j,k)+S(i,j,k+1)) & + + auxn(k) * ( Res (k+1) - Res (k) ) * (Res (k)+Res (k+1)) & + ) / dt +#endif + + end do + do k=1,kmax + phymix_S(i,j,k) = phymix_S(i,j,k) & + - ( hn(i,j,k)*(Res(k)*Res(k) - S(i,j,k)*S(i,j,k))/dt & + + (fluxw(k) - fluxw(k-1) ) ) + end do + end if do k=1,kmax S(i,j,k)=Res(k) @@ -520,8 +765,55 @@ !$OMP END PARALLEL + if (allocated(counts_s).or.allocated(flags_s).or.allocated(h_s).or.allocated(hS_s).or.allocated(hS2_s).or.allocated(hpmS_s).or.allocated(hnmS_s).or.allocated(phymix_S_fwf_s).or.allocated(phymix_S_riv_s)) then + call tic(TIM_TEF) + if (allocated(counts_s)) counts_s = 0 + if (allocated(flags_s )) flags_s = 0 + if (allocated(h_s )) h_s = _ZERO_ + if (allocated(hS_s )) hS_s = _ZERO_ + if (allocated(hS2_s )) hS2_s = _ZERO_ + if (allocated(hpmS_s )) hpmS_s = _ZERO_ + if (allocated(hnmS_s )) hnmS_s = _ZERO_ + if (allocated(phymix_S_fwf_s)) phymix_S_fwf_s = _ZERO_ + if (allocated(phymix_S_riv_s)) phymix_S_riv_s = _ZERO_ + do j = jmin,jmax + do i = imin,imax + if (az(i,j) .eq. 1) then + call rebin(kmax,S(i,j,:),saltbins,Si_s,binidx) + binidx_kij(:,i,j) = binidx(:) + do k=1,kmax + l = binidx(k) + if (l .lt. 0) cycle + if (allocated(counts_s)) counts_s(i,j,l) = counts_s(i,j,l) + 1 + if (allocated(flags_s )) flags_s (i,j,l) = 1 + if (allocated(h_s )) h_s (i,j,l) = h_s (i,j,l) + hn(i,j,k) + if (allocated(hS_s )) hS_s (i,j,l) = hS_s (i,j,l) + hn(i,j,k)*S(i,j,k) + if (allocated(hS2_s )) hS2_s (i,j,l) = hS2_s (i,j,l) + hn(i,j,k)*S(i,j,k)*S(i,j,k) + if (allocated(hpmS_s)) hpmS_s(i,j,l) = hpmS_s(i,j,l) + phymix_S(i,j,k) + if (allocated(hnmS_s)) hnmS_s(i,j,l) = hnmS_s(i,j,l) + nummix_S(i,j,k) + end do + l = binidx(kmax) + if (l .ge. 0) then + if (allocated(phymix_S_fwf_s)) phymix_S_fwf_s(i,j,l) = phymix_S_fwf(i,j) + end if + end if + end do + end do + if (allocated(phymix_S_riv_s)) then + do n=1,nriverl + j = rj(n) + i = ri(n) + do k=1,kmax + l = binidx_kij(k,i,j) + if (l .lt. 0) cycle + phymix_S_riv_s(i,j,l) = phymix_S_riv_s(i,j,l) + rivmix(k,n) + end do + end do + end if + call toc(TIM_TEF) + end if - if (salt_check .ne. 0 .and. mod(n,abs(salt_check)) .eq. 0) then + if (salt_check .ne. 0 .and. mod(loop,abs(salt_check)) .eq. 0) then call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, & S,min_salt,max_salt,status) if (status .gt. 0) then diff --git a/src/3d/tracer_diffusion.F90 b/src/3d/tracer_diffusion.F90 index d97c8556..df19e5db 100644 --- a/src/3d/tracer_diffusion.F90 +++ b/src/3d/tracer_diffusion.F90 @@ -23,6 +23,7 @@ use variables_les, only: AmU_3d,AmV_3d use m3d, only: AH_min_method use variables_3d, only: AH_min + use advection_3d, only: ffluxu2_3d,ffluxv2_3d use getm_timers, only: tic,toc,TIM_TRACEDIFF !$ use omp_lib IMPLICIT NONE @@ -130,6 +131,8 @@ calc_phymix = .false. end if + ffluxu2_3d = _ZERO_ ; ffluxv2_3d = _ZERO_ + ! Note (KK): diffusion only within new layer heigts !$OMP PARALLEL DEFAULT(SHARED) & @@ -420,6 +423,8 @@ end do #ifndef SLICE_MODEL end do + ffluxu2_3d(:,:,k) = ffluxu2 + ffluxv2_3d(:,:,k) = ffluxv2 #endif end if diff --git a/src/3d/variables_3d.F90 b/src/3d/variables_3d.F90 index baefa555..4849cbec 100644 --- a/src/3d/variables_3d.F90 +++ b/src/3d/variables_3d.F90 @@ -111,13 +111,13 @@ ! {\tt init\_variables\_3d}) and cleanup (see {\tt clean\_variables\_3d}). ! ! !USES: + use parameters, only: rk use domain, only: imin,imax,jmin,jmax,kmax,az,bottfric_method,rdrag use waves , only: waveforcing_method,waves_method,NO_WAVES,WAVES_VF use waves , only: waves_bbl_method,NO_WBBL IMPLICIT NONE ! ! !PUBLIC DATA MEMBERS: - integer, parameter :: rk = kind(_ONE_) REALTYPE :: dt,cnpar=0.9 REALTYPE :: avmback=_ZERO_,avhback=_ZERO_ logical :: deformC_3d=.false. @@ -129,6 +129,8 @@ logical :: save_Sfluxu=.false. logical :: save_Sfluxv=.false. logical :: save_Sfluxw=.false. + logical :: save_Sfluxu2=.false. + logical :: save_Sfluxv2=.false. logical :: save_Tfluxu=.false. logical :: save_Tfluxv=.false. logical :: save_Tfluxw=.false. @@ -136,8 +138,12 @@ logical :: save_phydis_3d=.false. logical :: save_nummix_S=.false. logical :: save_phymix_S=.false. + logical :: save_phymix_S_fwf=.false. logical :: save_nummix_T=.false. logical :: save_phymix_T=.false. + integer :: saltbins = 0 + REALTYPE :: saltbinmin=_ZERO_,saltbinmax=40.0_rk + integer :: id_dim_salt ! #ifdef STATIC #include "static_3d.h" @@ -160,6 +166,15 @@ REALTYPE,dimension(:,:,:),pointer,contiguous :: Tfluxu=>null() REALTYPE,dimension(:,:,:),pointer,contiguous :: Tfluxv=>null() REALTYPE,dimension(:,:,:),pointer,contiguous :: Tfluxw=>null() + REALTYPE,dimension(:,:,:),allocatable :: Sfluxu2 + REALTYPE,dimension(:,:,:),allocatable :: Sfluxv2 + + REALTYPE,dimension(:) ,allocatable :: Si_s + REALTYPE,dimension(:,:,:),allocatable :: counts_s,flags_s,h_s,hS_s,hS2_s,hpmS_s,hnmS_s + REALTYPE,dimension(:,:,:),allocatable :: hu_s,uu_s,Sfluxu_s,S2fluxu_s + REALTYPE,dimension(:,:,:),allocatable :: hv_s,vv_s,Sfluxv_s,S2fluxv_s + REALTYPE,dimension(:,:,:),allocatable :: fwf_s,fwfS2_s,phymix_S_fwf_s + REALTYPE,dimension(:,:,:),allocatable :: rvol_s,phymix_S_riv_s ! KK-TODO: make *dis_3d allocatable REALTYPE,dimension(:,:,:),pointer,contiguous :: numdis_3d=>null() @@ -168,6 +183,7 @@ REALTYPE,dimension(:,:,:),pointer,contiguous :: nummix_T=>null() REALTYPE,dimension(:,:,:),pointer,contiguous :: phymix_S=>null() REALTYPE,dimension(:,:,:),pointer,contiguous :: phymix_T=>null() + REALTYPE, dimension(:,:), allocatable :: phymix_S_fwf REALTYPE,public,dimension(:,:),allocatable,target :: sst,sss @@ -206,6 +222,27 @@ integer :: size3d_field integer :: mem3d integer :: preadapt + + logical, private :: counts_s_used=.false. + logical, private :: flags_s_used=.false. + logical :: h_s_used=.false. + logical :: hS_s_used=.false. + logical :: hS2_s_used=.false. + logical :: hpmS_s_used=.false. + logical :: hnmS_s_used=.false. + logical, private :: hu_s_used=.false. + logical, private :: hv_s_used=.false. + logical :: uu_s_used=.false. + logical :: vv_s_used=.false. + logical :: Sfluxu_s_used=.false. + logical :: Sfluxv_s_used=.false. + logical, private :: S2fluxu_s_used=.false. + logical, private :: S2fluxv_s_used=.false. + logical :: fwf_s_used=.false. + logical, private :: fwfS2_s_used=.false. + logical :: phymix_S_fwf_s_used=.false. + logical :: rvol_s_used=.false. + logical :: phymix_S_riv_s_used=.false. ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Hans Burchard @@ -390,6 +427,11 @@ ! !DESCRIPTION: ! ! !LOCAL VARIABLES: + logical :: calc_phymix_S=.false. + logical :: calc_nummix_S=.false. + logical :: calc_phymix_S_fwf=.false. + logical :: calc_Sfluxu,calc_Sfluxv + logical :: calc_Sfluxu2,calc_Sfluxv2 integer :: rc ! !EOP @@ -401,6 +443,11 @@ write(debug,*) 'postinit_variables_3d() # ',Ncall #endif + calc_Sfluxu = save_Sfluxu + calc_Sfluxv = save_Sfluxv + calc_Sfluxu2 = save_Sfluxu2 + calc_Sfluxv2 = save_Sfluxv2 + ! must be in postinit because flags are set init_getm_fabm if (deformC_3d) then allocate(dudxC_3d(I3DFIELD),stat=rc) @@ -463,12 +510,129 @@ #endif end if - if (save_Sfluxu) then + !if (update_salt .and. saltbins.gt.0) then + if (counts_s_used) then + allocate(counts_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (counts_s)' + counts_s = 0 + end if + if (flags_s_used) then + allocate(flags_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (flags_s)' + flags_s = 0 + end if + if (h_s_used) then + allocate(h_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (h_s)' + h_s = _ZERO_ + end if + if (hS_s_used) then + allocate(hS_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hS_s)' + hS_s = _ZERO_ + end if + if (hS2_s_used) then + allocate(hS2_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hS2_s)' + hS2_s = _ZERO_ + end if + if (hpmS_s_used) then + calc_phymix_S = .true. + allocate(hpmS_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hpmS_s)' + hpmS_s = _ZERO_ + end if + if (hnmS_s_used) then + calc_nummix_S = .true. + allocate(hnmS_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hnmS_s)' + hnmS_s = _ZERO_ + end if + if (hu_s_used) then + calc_Sfluxu = .true. + allocate(hu_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hu_s)' + hu_s = _ZERO_ + end if + if (hv_s_used) then + calc_Sfluxv = .true. + allocate(hv_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (hv_s)' + hv_s = _ZERO_ + end if + if (uu_s_used) then + calc_Sfluxu = .true. + allocate(uu_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (uu_s)' + uu_s = _ZERO_ + end if + if (vv_s_used) then + calc_Sfluxv = .true. + allocate(vv_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (vv_s)' + vv_s = _ZERO_ + end if + if (Sfluxu_s_used) then + calc_Sfluxu = .true. + allocate(Sfluxu_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxu_s)' + Sfluxu_s = _ZERO_ + end if + if (S2fluxu_s_used) then + calc_Sfluxu = .true. + calc_Sfluxu2 = .true. + allocate(S2fluxu_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (S2fluxu_s)' + S2fluxu_s = _ZERO_ + end if + if (Sfluxv_s_used) then + calc_Sfluxv = .true. + allocate(Sfluxv_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxv_s)' + Sfluxv_s = _ZERO_ + end if + if (S2fluxv_s_used) then + calc_Sfluxv = .true. + calc_Sfluxv2 = .true. + allocate(S2fluxv_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (S2fluxv_s)' + S2fluxv_s = _ZERO_ + end if + if (fwf_s_used) then + allocate(fwf_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (fwf_s)' + fwf_s = _ZERO_ + end if + if (fwfS2_s_used) then + allocate(fwfS2_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (fwfS2_s)' + fwfS2_s = _ZERO_ + end if + if (phymix_S_fwf_s_used) then + calc_phymix_S_fwf = .true. + allocate(phymix_S_fwf_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_fwf_s)' + phymix_S_fwf_s = _ZERO_ + end if + if (rvol_s_used) then + allocate(rvol_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (rvol_s)' + rvol_s = _ZERO_ + end if + if (phymix_S_riv_s_used) then + calc_phymix_S = .true. + allocate(phymix_S_riv_s(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_riv_s)' + phymix_S_riv_s = _ZERO_ + end if + !end if + + if (calc_Sfluxu) then allocate(Sfluxu(I3DFIELD),stat=rc) if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxu)' Sfluxu = _ZERO_ end if - if (save_Sfluxv) then + if (calc_Sfluxv) then allocate(Sfluxv(I3DFIELD),stat=rc) if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxv)' Sfluxv = _ZERO_ @@ -478,6 +642,17 @@ if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxw)' Sfluxw = _ZERO_ end if + if (calc_Sfluxu2) then + allocate(Sfluxu2(I3DFIELD),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxu2)' + Sfluxu2 = _ZERO_ + end if + if (calc_Sfluxv2) then + allocate(Sfluxv2(I3DFIELD),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (Sfluxv2)' + Sfluxv2 = _ZERO_ + end if + if (save_Tfluxu) then allocate(Tfluxu(I3DFIELD),stat=rc) if (rc /= 0) stop 'postinit_3d: Error allocating memory (Tfluxu)' @@ -527,18 +702,24 @@ nummix_T = _ZERO_ end if - if (save_phymix_S) then + if (save_phymix_S .or. calc_phymix_S) then save_phymix_3d=.true. allocate(phymix_S(I3DFIELD),stat=rc) if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S)' phymix_S = _ZERO_ end if - if (save_nummix_S) then + if (save_nummix_S .or. calc_nummix_S) then allocate(nummix_S(I3DFIELD),stat=rc) if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S)' nummix_S = _ZERO_ end if + if (save_phymix_S_fwf .or. calc_phymix_S_fwf) then + allocate(phymix_S_fwf(I2DFIELD),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_fwf)' + phymix_S_fwf = _ZERO_ + end if + #ifdef DEBUG write(debug,*) 'Leaving postinit_variables_3d()' write(debug,*) @@ -561,6 +742,7 @@ use field_manager use meteo , only: fwf_method use variables_2d, only: fwf_int + use domain , only: nriverl IMPLICIT NONE ! ! !INPUT PARAMETERS: @@ -674,6 +856,8 @@ ! category - 3d if (runtype .ge. 2) then + call fm%register('Uavg', 'm2/s', 'transport in local x-direction (2D)', standard_name='', data2d=Uavg(_2D_FO_), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, part_of_state=.true.) + call fm%register('Vavg', 'm2/s', 'transport in local y-direction (2D)', standard_name='', data2d=Vavg(_2D_FO_), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, part_of_state=.true.) call fm%register('zc', 'm', 'center coordinate', standard_name='', dimensions=(/id_dim_z/),data3d=zc(_3D_FO_), update_interval=update_interval, category='grid', part_of_state=.false.) call fm%register('zcn', 'm', 'z', standard_name='', dimensions=(/id_dim_z/), data3d=zcn(_3D_FO_), update_interval=update_interval, category='grid') call fm%register('hn', 'm', 'layer thickness', standard_name='cell_thickness', dimensions=(/id_dim_z/),data3d=hn(_3D_FO_), update_interval=update_interval, category='grid', part_of_state=.true.) @@ -737,7 +921,41 @@ call fm%register('Sfluxu', 'g/kg*m3/s', 'salt flux in local x-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Sfluxu) call fm%register('Sfluxv', 'g/kg*m3/s', 'salt flux in local y-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Sfluxv) call fm%register('Sfluxw', 'g/kg*m/s', 'vertical salt flux', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Sfluxw) + call fm%register('Sfluxu2', '(g/kg)**2*m3/s', 'salt variance flux in local x-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Sfluxu2) + call fm%register('Sfluxv2', '(g/kg)**2*m3/s', 'salt variance flux in local y-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Sfluxv2) + + if (saltbins .gt. 0) then + !call fm%register_dimension('salt_s', saltbins+1, global_length=saltbins, offset=-1, id=id_dim_salt) + call fm%register_dimension('salt_s', saltbins+1, newid=id_dim_salt) + call fm%register('salt_s', 'g/kg', 'interface salinity (bin)', dimensions=(/id_dim_salt/), no_default_dimensions=.true., coordinate_dimension=id_dim_salt, category='TEF', output_level=output_level_debug) + call fm%register('counts_s', '1', 'bin counts', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', output_level=output_level_debug, used=counts_s_used) + call fm%register('flags_s', '1', 'bin flags', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', output_level=output_level_debug, used=flags_s_used) + call fm%register('h_s', 'm', 'bin thickness', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', output_level=output_level_debug, used=h_s_used) + call fm%register('hS_s', 'm*g/kg', 'salt content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', output_level=output_level_debug, used=hS_s_used) + call fm%register('hS2_s', 'm*(g/kg)**2', 'salt2 content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', output_level=output_level_debug, used=hS2_s_used) + call fm%register('hpmS_s', 'm*(g/kg)**2/s', 'phymixS content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=hpmS_s_used) + call fm%register('hnmS_s', 'm*(g/kg)**2/s', 'nummixS content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=hnmS_s_used) + call fm%register('hu_s', 'm' , 'bin thickness at U-points', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=hu_s_used) + call fm%register('hv_s', 'm' , 'bin thickness at V-points', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=hv_s_used) + call fm%register('uu_s', 'm2/s', 'transport in local x-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=uu_s_used) + call fm%register('vv_s', 'm2/s', 'transport in local y-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=vv_s_used) + call fm%register('Sfluxu_s', 'g/kg*m3/s', 'salt flux in local x-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=Sfluxu_s_used) + call fm%register('Sfluxv_s', 'g/kg*m3/s', 'salt flux in local y-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=Sfluxv_s_used) + call fm%register('S2fluxu_s', '(g/kg)**2*m3/s', 'salt squared flux in local x-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=S2fluxu_s_used) + call fm%register('S2fluxv_s', '(g/kg)**2*m3/s', 'salt squared flux in local y-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=S2fluxv_s_used) + if (fwf_method .ge. 1) then + call fm%register('fwf_s', 'm', 'fresh water fluxes (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=fwf_s_used) + call fm%register('fwfS2_s', 'm*(g/kg)**2', 'fwfS2_s', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=fwfS2_s_used) + call fm%register('phymix_S_fwf_s', 'm*(g/kg)**2/s', 'physical mixing content of salinity (fwf) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=phymix_S_fwf_s_used) + end if + if (nriverl .gt. 0) then + call fm%register('rvol_s', 'm**3', 'river discharge (3D) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=rvol_s_used) + call fm%register('phymix_S_riv_s', 'm*(g/kg)**2/s', 'physical mixing content of salinity (riv) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEF', add_to_mean_before_save=.true., output_level=output_level_debug, used=phymix_S_riv_s_used) + end if + + end if end if + if (update_temp) then call fm%register('Tfluxu', 'degree_C*m3/s', 'temp flux in local x-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Tfluxu) call fm%register('Tfluxv', 'degree_C*m3/s', 'temp flux in local y-direction', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_Tfluxv) @@ -751,6 +969,9 @@ if (update_salt) then call fm%register('nummix_S', 'm*(g/kg)**2/s', 'numerical mixing content of salinity', standard_name='', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_nummix_S) call fm%register('phymix_S', 'm*(g/kg)**2/s', 'physical mixing content of salinity' , standard_name='', dimensions=(/id_dim_z/), update_interval=update_interval, category='3d', add_to_mean_before_save=.true., output_level=output_level_debug, used=save_phymix_S) + if (fwf_method .ge. 1) then + call fm%register('phymix_S_fwf', 'm*(g/kg)**2/s', 'physical mixing content of salinity (fwf)', standard_name='', update_interval=update_interval, category='3d', fill_value=-9999.0_rk, add_to_mean_before_save=.true., output_level=output_level_debug, used=save_phymix_S_fwf) + end if end if return @@ -787,6 +1008,31 @@ if (associated(Sfluxu)) call fm%send_data('Sfluxu', Sfluxu(_3D_FO_)) if (associated(Sfluxv)) call fm%send_data('Sfluxv', Sfluxv(_3D_FO_)) if (associated(Sfluxw)) call fm%send_data('Sfluxw', Sfluxw(_3D_FO_)) + if (allocated(Sfluxu2)) call fm%send_data('Sfluxu2', Sfluxu2(_3D_FO_)) + if (allocated(Sfluxv2)) call fm%send_data('Sfluxv2', Sfluxv2(_3D_FO_)) + + if (allocated(Si_s )) call fm%send_data('salt_s', Si_s) + if (allocated(counts_s)) call fm%send_data('counts_s', counts_s(_2D_FO_,:)) + if (allocated(flags_s )) call fm%send_data('flags_s' , flags_s (_2D_FO_,:)) + if (allocated(h_s )) call fm%send_data('h_s' , h_s (_2D_FO_,:)) + if (allocated(hS_s )) call fm%send_data('hS_s' , hS_s (_2D_FO_,:)) + if (allocated(hS2_s )) call fm%send_data('hS2_s' , hS2_s (_2D_FO_,:)) + if (allocated(hpmS_s)) call fm%send_data('hpmS_s', hpmS_s(_2D_FO_,:)) + if (allocated(hnmS_s)) call fm%send_data('hnmS_s', hnmS_s(_2D_FO_,:)) + if (allocated(hu_s )) call fm%send_data('hu_s' , hu_s (_2D_FO_,:)) + if (allocated(hv_s )) call fm%send_data('hv_s' , hv_s (_2D_FO_,:)) + if (allocated(uu_s )) call fm%send_data('uu_s' , uu_s (_2D_FO_,:)) + if (allocated(vv_s )) call fm%send_data('vv_s' , vv_s (_2D_FO_,:)) + if (allocated(Sfluxu_s)) call fm%send_data('Sfluxu_s', Sfluxu_s(_2D_FO_,:)) + if (allocated(Sfluxv_s)) call fm%send_data('Sfluxv_s', Sfluxv_s(_2D_FO_,:)) + if (allocated(S2fluxu_s)) call fm%send_data('S2fluxu_s', S2fluxu_s(_2D_FO_,:)) + if (allocated(S2fluxv_s)) call fm%send_data('S2fluxv_s', S2fluxv_s(_2D_FO_,:)) + if (allocated(fwf_s )) call fm%send_data('fwf_s' , fwf_s (_2D_FO_,:)) + if (allocated(fwfS2_s)) call fm%send_data('fwfS2_s', fwfS2_s(_2D_FO_,:)) + if (allocated(phymix_S_fwf_s)) call fm%send_data('phymix_S_fwf_s', phymix_S_fwf_s(_2D_FO_,:)) + if (allocated(rvol_s)) call fm%send_data('rvol_s', rvol_s(_2D_FO_,:)) + if (allocated(phymix_S_riv_s)) call fm%send_data('phymix_S_riv_s', phymix_S_riv_s(_2D_FO_,:)) + if (associated(Tfluxu)) call fm%send_data('Tfluxu', Tfluxu(_3D_FO_)) if (associated(Tfluxv)) call fm%send_data('Tfluxv', Tfluxv(_3D_FO_)) if (associated(Tfluxw)) call fm%send_data('Tfluxw', Tfluxw(_3D_FO_)) @@ -796,6 +1042,7 @@ if (associated(phymix_T )) call fm%send_data('phymix_T' , phymix_T (_3D_FO_)) if (associated(nummix_S )) call fm%send_data('nummix_S' , nummix_S (_3D_FO_)) if (associated(phymix_S )) call fm%send_data('phymix_S' , phymix_S (_3D_FO_)) + if (allocated(phymix_S_fwf)) call fm%send_data('phymix_S_fwf', phymix_S_fwf(_2D_FO_)) return end subroutine finalize_register_3d_variables diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a05c4a2f..f923c2aa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -255,6 +255,7 @@ add_library(futils OBJECT futils/parameters.F90 futils/pos.F90 futils/read_par_setup.F90 + futils/rebin.F90 futils/strip_string.F90 futils/time.F90 futils/to_2d_u.F90 diff --git a/src/domain/domain.F90 b/src/domain/domain.F90 index 50f1cfc6..38626ebd 100644 --- a/src/domain/domain.F90 +++ b/src/domain/domain.F90 @@ -95,6 +95,9 @@ logical :: need_2d_bdy_v = .false. logical :: need_3d_bdy = .false. + integer :: nriverl=0 ! local number of river cells + integer, dimension(:), allocatable :: ri,rj + REALTYPE :: cori= _ZERO_ integer :: bottfric_method=2 diff --git a/src/futils/Makefile b/src/futils/Makefile index 00722e2d..72b29841 100644 --- a/src/futils/Makefile +++ b/src/futils/Makefile @@ -61,6 +61,7 @@ ${LIB}(to_v.o) \ ${LIB}(to_w.o) \ ${LIB}(c2x.o) \ ${LIB}(check_3d_fields.o) \ +$(LIB)(rebin.o) \ ${LIB}(filter_1d.o) all: modules objects diff --git a/src/futils/getm_timers.F90 b/src/futils/getm_timers.F90 index 0be2f140..2296cfd5 100644 --- a/src/futils/getm_timers.F90 +++ b/src/futils/getm_timers.F90 @@ -87,6 +87,7 @@ integer, parameter :: TIM_ADV3DH = 103 ! 3d advection halo parts integer, parameter :: TIM_CHECK3DF = 104 ! check_3d_fields integer, parameter :: TIM_MIXANALYSIS = 105 ! (numerical) mixing analysis + integer, parameter :: TIM_TEF = 108 ! TEF analysis integer, parameter :: TIM_HALO2D = 110 ! do halo 2d (initialize comm) integer, parameter :: TIM_HALO3D = 111 ! do halo 3d (initialize comm) integer, parameter :: TIM_HALOWAIT = 112 ! wait_halo (2d+3d both) @@ -241,6 +242,7 @@ timernames(TIM_ADV3D) = ' sum do_advection_3d' timernames(TIM_INTEGR3D) = 'integrate_3d other' timernames(TIM_MIXANALYSIS) = 'numerical mixing analysis' + timernames(TIM_TEF) = ' TEF analysis' ! We only really want to display halo-stuff if we compile for parallel: #ifdef GETM_PARALLEL diff --git a/src/futils/rebin.F90 b/src/futils/rebin.F90 new file mode 100644 index 00000000..672cb2c5 --- /dev/null +++ b/src/futils/rebin.F90 @@ -0,0 +1,71 @@ +#include "cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: rebin - +! \label{sec-rebin} +! +! !INTERFACE: + subroutine rebin(nlev,flev,nbins,fibins,binidx) +! +! !DESCRIPTION: +! +! !USES: + IMPLICIT NONE +! +! !INPUT PARAMETERS: + integer, intent(in) :: nlev,nbins + REALTYPE, dimension(0:nlev), intent(in) :: flev + REALTYPE, dimension(0:nbins), intent(in) :: fibins +! +! !INPUT/OUTPUT PARAMETERS: +! +! !OUTPUT PARAMETERS: + integer, dimension(0:nlev), intent(out) :: binidx +! +! !REVISION HISTORY: +! Original author(s): Knut Klingbeil +! +! !LOCAL VARIABLES: + integer :: k,l,l0 +!EOP +!----------------------------------------------------------------------- +!BOC +#ifdef DEBUG + integer, save :: Ncall = 0 + Ncall = Ncall+1 + write(debug,*) 'rebin() # ',Ncall +#endif + + l0 = nbins + do k=1,nlev +! fibins strictly monotonically increasing with l + if (flev(k) .gt. fibins(nbins)) then + binidx(k) = -1 + cycle + end if +#if 1 +! equidistant bins + l = max( 0 , int( ceiling( (flev(k)-fibins(0))/(fibins(nbins)-fibins(0))*nbins ) ) ) +#else + if (flev(k) .gt. fibins( l0 )) l0 = nbins ! reset to start from highest bin again + do l = l0,1,-1 + if (flev(k) .gt. fibins(l-1)) then + l0 = l + exit + end if + end do +#endif + binidx(k) = l ! here we also map to index 0 + end do + +#ifdef DEBUG + write(debug,*) 'Leaving rebin()' + write(debug,*) +#endif + return + end subroutine rebin +!EOC +!----------------------------------------------------------------------- +! Copyright (C) 2019 - Knut Klingbeil ! +!----------------------------------------------------------------------- diff --git a/src/getm/compilation_options.F90 b/src/getm/compilation_options.F90 index ddfd0547..5e975472 100644 --- a/src/getm/compilation_options.F90 +++ b/src/getm/compilation_options.F90 @@ -166,6 +166,9 @@ #ifdef SONG_WRIGHT LEVEL1 'SONG_WRIGHT' #endif +#ifdef _MOMENTUM_TERMS_ + LEVEL1 '_MOMENTUM_TERMS_' +#endif #ifdef NO_TIMERS LEVEL1 'NO_TIMERS' #endif diff --git a/src/getm/initialise.F90 b/src/getm/initialise.F90 index 1c42a7e5..70fc7c01 100644 --- a/src/getm/initialise.F90 +++ b/src/getm/initialise.F90 @@ -371,8 +371,6 @@ call init_les(runtype) - call init_output_processing() - call init_register_all_variables(runtype) #ifdef _NEW_GOTM_ @@ -393,7 +391,12 @@ call init_output(runid,title,start,runtype,dryrun,myid,MinN,MaxN,save_initial) call do_register_all_variables(runtype) - if (list_variables) call fm%list() + if (list_variables) then + call fm%list() + stop 0 + end if + + call init_output_processing() close(NAMLST) @@ -524,7 +527,7 @@ #else call output_manager_prepare_save(julianday, int(fsecondsofday), microseconds, int(MinN-1,kind=timestepkind)) #endif - call do_output_processing() + call do_output_processing(MinN-1) #ifdef _NEW_GOTM_ call output_manager_save(julianday, int(fsecondsofday), microseconds, MinN-1)) #else diff --git a/src/getm/integration.F90 b/src/getm/integration.F90 index 76f18cea..0a9ee793 100644 --- a/src/getm/integration.F90 +++ b/src/getm/integration.F90 @@ -98,10 +98,11 @@ use variables_3d, only: sseo,ssen,ho,hn #ifndef NO_BAROCLINIC use variables_3d, only: rho,T,S + use variables_3d, only: phymix_T,phymix_S #endif use rivers, only: do_rivers #ifdef _FABM_ - use getm_fabm, only: fabm_calc,do_getm_fabm + use getm_fabm, only: fabm_calc,do_getm_fabm, pa_phymix_fabm_pel #endif #ifdef GETM_BIO use bio, only: bio_calc @@ -128,6 +129,7 @@ ! ! !LOCAL VARIABLES logical :: do_3d=.false. + integer :: i character(8) :: d_ character(10) :: t_ !EOP @@ -182,6 +184,16 @@ if (do_3d) then sseo = ssen ! true sseo (without rivers and fwf) ho = hn ! true ho (without rivers and fwf) +! phymix is only reset if rivers are present + if (associated(phymix_T)) phymix_T = _ZERO_ + if (associated(phymix_S)) phymix_S = _ZERO_ +#ifdef _FABM_ + if (allocated(pa_phymix_fabm_pel)) then + do i=1,size(pa_phymix_fabm_pel) + if (associated(pa_phymix_fabm_pel(i)%p3d)) pa_phymix_fabm_pel(i)%p3d = _ZERO_ + end do + end if +#endif end if call do_rivers(n,do_3d) if (do_3d) then @@ -216,7 +228,7 @@ #endif call toc(TIM_FLEX_OUTPUT) call tic(TIM_OUTPUT_PROC) - call do_output_processing() + call do_output_processing(n) call toc(TIM_OUTPUT_PROC) call tic(TIM_FLEX_OUTPUT) #ifdef _NEW_GOTM_ diff --git a/src/getm/register_all_variables.F90 b/src/getm/register_all_variables.F90 index f861e943..8fff7f11 100644 --- a/src/getm/register_all_variables.F90 +++ b/src/getm/register_all_variables.F90 @@ -197,7 +197,7 @@ #if 0 call register_diagnostic_variables() #endif - call register_processed_variables(fm,runtype) + call register_processed_variables(fm,runtype,M) return end subroutine do_register_all_variables @@ -296,8 +296,8 @@ call fm%register('yx', 'm', 'y at X-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=yx(_2D_FO_), category='domain', output_level=output_level_debug) end if if (grid_type .eq. 4) then - call fm%register('lonx', trim(lonunits), trim(lonlongname)//' at X-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=lonx(_2D_W_), category='domain', output_level=output_level_debug) - call fm%register('latx', trim(latunits), trim(latlongname)//' at X-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=latx(_2D_W_), category='domain', output_level=output_level_debug) + call fm%register('lonx', trim(lonunits), trim(lonlongname)//' at X-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=lonx(_2D_FO_), category='domain', output_level=output_level_debug) + call fm%register('latx', trim(latunits), trim(latlongname)//' at X-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=latx(_2D_FO_), category='domain', output_level=output_level_debug) end if call fm%register('bathymetry', 'm', 'bathymetry', standard_name='bathymetry', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., fill_value=-10._rk, data2d=H(_2D_FO_), category='domain',output_level=output_level_required) @@ -310,11 +310,11 @@ call fm%register('areaC', 'm2', 'area at T-points', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=areaC(_2D_FO_), category="metrics", output_level=output_level_debug) if (grid_type .ne. 2) then - call fm%register('convc', 'degree', 'clockwise angle between geographical and local axes', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=convc(_2D_W_), category="domain", output_level=output_level_debug) + call fm%register('convc', 'degree', 'clockwise angle between geographical and local axes', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=convc(_2D_FO_), category="domain", output_level=output_level_debug) end if if (grid_type.eq.3 .or. grid_type.eq.4) then - call fm%register('gtlc', 'degree', 'clockwise angle between global and local axes', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=gtlc(_2D_W_), category="domain", output_level=output_level_debug) + call fm%register('gtlc', 'degree', 'clockwise angle between global and local axes', dimensions=(/id_dim_lon,id_dim_lat/), no_default_dimensions=.true., data2d=gtlc(_2D_FO_), category="domain", output_level=output_level_debug) end if return diff --git a/src/ncdf/ncdf_rivers.F90 b/src/ncdf/ncdf_rivers.F90 index f31e9b0d..5ff65d34 100644 --- a/src/ncdf/ncdf_rivers.F90 +++ b/src/ncdf/ncdf_rivers.F90 @@ -298,8 +298,8 @@ if (err .ne. NF90_NOERR) go to 10 do m=1,river_split(ni) river_flow(ni+m-1) = river_factor*x(1) - river_salt(ni+m-1) = salt_missing - river_temp(ni+m-1) = temp_missing + !river_salt(ni+m-1) = salt_missing + !river_temp(ni+m-1) = temp_missing end do if ( r_salt(nn) .eq. 1 ) then err = nf90_get_var(ncid,salt_id(nn),x,start,edges) diff --git a/src/output/output_processing.F90 b/src/output/output_processing.F90 index fd25039e..41177774 100644 --- a/src/output/output_processing.F90 +++ b/src/output/output_processing.F90 @@ -11,6 +11,7 @@ ! This modules serves as a container for processing output variables. ! !USES: + use parameters, only: rk use domain, only: imin,imax,jmin,jmax,kmax use field_manager IMPLICIT NONE @@ -29,13 +30,53 @@ REALTYPE, dimension(:,:), allocatable, target :: u2d_destag, v2d_destag REALTYPE, dimension(:,:,:), allocatable, target :: u3d_destag, v3d_destag +#ifndef NO_BAROCLINIC + REALTYPE,dimension(:,:,:),allocatable,target :: h_s_mean, hS_s_mean, hS2_s_mean + REALTYPE,dimension(:,:,:),allocatable,target :: uu_s_mean, vv_s_mean + REALTYPE,dimension(:,:,:),allocatable,target :: Sfluxu_s_mean, Sfluxv_s_mean + REALTYPE,dimension(:,:,:),allocatable,target :: hpmS_s_mean, hnmS_s_mean + REALTYPE,dimension(:,:,:),allocatable,target :: fwf_s_int, phymix_S_fwf_s_mean + REALTYPE,dimension(:,:,:),allocatable,target :: rvol_s_int, phymix_S_riv_s_mean +#endif + logical :: u2d_used, v2d_used logical :: u3d_used=.false., v3d_used=.false. + logical :: u2d_destag_used, v2d_destag_used logical, target :: u2d_now, v2d_now logical, target :: u3d_now=.false., v3d_now=.false. - logical, target:: u2d_destag_use, v2d_destag_use + logical, target :: u2d_destag_now, v2d_destag_now logical, target:: u3d_destag_use, v3d_destag_use - integer, parameter :: rk = kind(_ONE_) + +#ifndef NO_BAROCLINIC + logical :: process_TEFmean=.false. + logical :: h_s_mean_used=.false. + logical :: hS_s_mean_used=.false. + logical :: hS2_s_mean_used=.false. + logical :: uu_s_mean_used=.false. + logical :: vv_s_mean_used=.false. + logical :: Sfluxu_s_mean_used=.false. + logical :: Sfluxv_s_mean_used=.false. + logical :: hpmS_s_mean_used=.false. + logical :: hnmS_s_mean_used=.false. + logical :: fwf_s_int_used=.false. + logical :: phymix_S_fwf_s_mean_used=.false. + logical :: rvol_s_int_used=.false. + logical :: phymix_S_riv_s_mean_used=.false. + logical, target :: h_s_mean_now=.false. + logical, target :: hS_s_mean_now=.false. + logical, target :: hS2_s_mean_now=.false. + logical, target :: uu_s_mean_now=.false. + logical, target :: vv_s_mean_now=.false. + logical, target :: Sfluxu_s_mean_now=.false. + logical, target :: Sfluxv_s_mean_now=.false. + logical, target :: hpmS_s_mean_now=.false. + logical, target :: hnmS_s_mean_now=.false. + logical, target :: fwf_s_int_now=.false. + logical, target :: phymix_S_fwf_s_mean_now=.false. + logical, target :: rvol_s_int_now=.false. + logical, target :: phymix_S_riv_s_mean_now=.false. +#endif + ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Hans Burchard @@ -54,6 +95,15 @@ subroutine init_output_processing ! ! !USES: +#ifndef NO_BAROCLINIC + use variables_3d, only: saltbins + use variables_3d, only: h_s_used, hS_s_used, hS2_s_used + use variables_3d, only: uu_s_used, vv_s_used + use variables_3d, only: Sfluxu_s_used, Sfluxv_s_used + use variables_3d, only: hpmS_s_used, hnmS_s_used + use variables_3d, only: fwf_s_used, phymix_S_fwf_s_used + use variables_3d, only: rvol_s_used, phymix_S_riv_s_used +#endif IMPLICIT NONE ! ! !DESCRIPTION: @@ -70,19 +120,99 @@ !------------------------------------------------------------------------- !BOC - allocate(u2d_destag(E2DFIELD),stat=rc) - if (rc /= 0) stop 'init_output_processing: Error allocating memory (u2d_destag)' - u2d_destag = 0._rk - allocate(v2d_destag(E2DFIELD),stat=rc) - if (rc /= 0) stop 'init_output_processing: Error allocating memory (v2d_destag)' - v2d_destag = 0._rk - #if 0 allocate(u3d_destag(I3DFIELD),stat=rc) if (rc /= 0) stop 'init_output_processing: Error allocating memory (u3d_destag)' allocate(v3d_destag(I3DFIELD),stat=rc) if (rc /= 0) stop 'init_output_processing: Error allocating memory (v3d_destag)' #endif + +! Note (KK): this needs to be called +! AFTER register_processed_variables() +! and BEFORE postinit_variables_3d() !!! +#ifndef NO_BAROCLINIC + if (process_TEFmean) then + if (h_s_mean_used) then + h_s_used = .true. + allocate(h_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (h_s_mean)' + h_s_mean = _ZERO_ + end if + if (hS_s_mean_used) then + hS_s_used = .true. + allocate(hS_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (hS_s_mean)' + hS_s_mean = _ZERO_ + end if + if (hS2_s_mean_used) then + hS2_s_used = .true. + allocate(hS2_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (hS2_s_mean)' + hS2_s_mean = _ZERO_ + end if + if (uu_s_mean_used) then + uu_s_used = .true. + allocate(uu_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (uu_s_mean)' + uu_s_mean = _ZERO_ + end if + if (vv_s_mean_used) then + vv_s_used = .true. + allocate(vv_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (vv_s_mean)' + vv_s_mean = _ZERO_ + end if + if (Sfluxu_s_mean_used) then + Sfluxu_s_used = .true. + allocate(Sfluxu_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (Sfluxu_s_mean)' + Sfluxu_s_mean = _ZERO_ + end if + if (Sfluxv_s_mean_used) then + Sfluxv_s_used = .true. + allocate(Sfluxv_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (Sfluxv_s_mean)' + Sfluxv_s_mean = _ZERO_ + end if + if (hpmS_s_mean_used) then + hpmS_s_used = .true. + allocate(hpmS_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (hpmS_s_mean)' + hpmS_s_mean = _ZERO_ + end if + if (hnmS_s_mean_used) then + hnmS_s_used = .true. + allocate(hnmS_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (hnmS_s_mean)' + hnmS_s_mean = _ZERO_ + end if + if (fwf_s_int_used) then + fwf_s_used = .true. + allocate(fwf_s_int(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_variables_3d: Error allocating memory (fwf_s_int)' + fwf_s_int = _ZERO_ + end if + if (phymix_S_fwf_s_mean_used) then + phymix_S_fwf_s_used = .true. + allocate(phymix_S_fwf_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_fwf_s_mean)' + phymix_S_fwf_s_mean = _ZERO_ + end if + if (rvol_s_int_used) then + rvol_s_used = .true. + allocate(rvol_s_int(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (rvol_s_int)' + rvol_s_int = _ZERO_ + end if + if (phymix_S_riv_s_mean_used) then + phymix_S_riv_s_used = .true. + allocate(phymix_S_riv_s_mean(I2DFIELD,0:saltbins),stat=rc) + if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_riv_s_mean)' + phymix_S_riv_s_mean = _ZERO_ + end if + end if +#endif + return end subroutine init_output_processing !EOC @@ -93,26 +223,37 @@ ! !ROUTINE: register_processed_variables() ! ! !INTERFACE: - subroutine register_processed_variables(fm,runtype) + subroutine register_processed_variables(fm,runtype,M) ! ! !DESCRIPTION: ! ! !USES: +#ifndef NO_BAROCLINIC + use m3d, only: update_salt + use variables_3d, only: saltbins + use variables_3d, only: id_dim_salt +#endif + use meteo , only: fwf_method + use domain , only: nriverl IMPLICIT NONE ! ! !INPUT PARAMETERS: type (type_field_manager) :: fm integer, intent(in) :: runtype + integer, intent(in) :: M ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding & Jorn Bruggeman ! ! !LOCAL VARIABLES: + integer :: update_interval !EOP !----------------------------------------------------------------------- !BOC LEVEL2 'register_processed_variables()' + update_interval = M + call fm%register('u2d', 'm/s', 'velocity in local x-direction', standard_name='', fill_value=-9999._rk, category='velocities', output_level=output_level_debug, used=u2d_used, used_now=u2d_now) call fm%register('v2d', 'm/s', 'velocity in local y-direction', standard_name='', fill_value=-9999._rk, category='velocities', output_level=output_level_debug, used=v2d_used, used_now=v2d_now) @@ -124,8 +265,38 @@ #endif - call fm%register('u2d-destag', 'm/s', 'velocity in local x-direction(destag)', standard_name='', data2d=u2d_destag(_2D_FO_), fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used_now=u2d_destag_use) - call fm%register('v2d-destag', 'm/s', 'velocity in local y-direction(destag)', standard_name='', data2d=v2d_destag(_2D_FO_), fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used_now=v2d_destag_use) + call fm%register('u2d-destag', 'm/s', 'velocity in local x-direction(destag)', standard_name='', fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used=u2d_destag_used, used_now=u2d_destag_now) + call fm%register('v2d-destag', 'm/s', 'velocity in local y-direction(destag)', standard_name='', fill_value=-9999._rk, category='velocities',output_level=output_level_debug, used=v2d_destag_used, used_now=v2d_destag_now) + +#ifndef NO_BAROCLINIC + if (runtype .ge. 3) then + if (update_salt) then + if (saltbins .gt. 0) then + + process_TEFmean = .true. + + call fm%register('h_s_mean', 'm', 'mean bin thickness', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', output_level=output_level_debug, used=h_s_mean_used, used_now=h_s_mean_now) + call fm%register('hS_s_mean', 'm*g/kg', 'mean salt content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', output_level=output_level_debug, used=hS_s_mean_used, used_now=hS_s_mean_now) + call fm%register('hS2_s_mean', 'm*(g/kg)**2', 'mean salt2 content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', output_level=output_level_debug, used=hS2_s_mean_used, used_now=hS2_s_mean_now) + call fm%register('uu_s_mean', 'm2/s', 'mean transport in local x-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=uu_s_mean_used, used_now=uu_s_mean_now) + call fm%register('vv_s_mean', 'm2/s', 'mean transport in local y-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=vv_s_mean_used, used_now=vv_s_mean_now) + call fm%register('Sfluxu_s_mean', '(g/kg)*m3/s', 'mean salt flux in local x-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=Sfluxu_s_mean_used, used_now=Sfluxu_s_mean_now) + call fm%register('Sfluxv_s_mean', '(g/kg)*m3/s', 'mean salt flux in local y-direction (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=Sfluxv_s_mean_used, used_now=Sfluxv_s_mean_now) + call fm%register('hpmS_s_mean', 'm*(g/kg)**2/s', 'mean phymixS content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=hpmS_s_mean_used, used_now=hpmS_s_mean_now) + call fm%register('hnmS_s_mean', 'm*(g/kg)**2/s', 'mean nummixS content (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=hnmS_s_mean_used, used_now=hnmS_s_mean_now) + if (fwf_method .ge. 1) then + call fm%register('fwf_s_int', 'm', 'integrated fresh water fluxes (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=fwf_s_int_used, used_now=fwf_s_int_now) + call fm%register('phymix_S_fwf_s_mean', 'm*(g/kg)**2/s', 'mean physical mixing content of salinity (fwf) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=phymix_S_fwf_s_mean_used, used_now=phymix_S_fwf_s_mean_now) + end if + if (nriverl .gt. 0) then + call fm%register('rvol_s_int', 'm**3', 'integrated river discharge (3D) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=rvol_s_int_used, used_now=rvol_s_int_now) + call fm%register('phymix_S_riv_s_mean', 'm*(g/kg)**2/s', 'mean physical mixing content of salinity (riv) (bin)', dimensions=(/id_dim_salt/), update_interval=update_interval, category='TEFmean', add_to_mean_before_save=.true., output_level=output_level_debug, used=phymix_S_riv_s_mean_used, used_now=phymix_S_riv_s_mean_now) + end if + + end if + end if + end if +#endif return end subroutine register_processed_variables @@ -188,6 +359,39 @@ end if #endif + if (u2d_destag_used) then + allocate(u2d_destag(E2DFIELD),stat=rc) + if (rc /= 0) stop 'finalize_register_processed_variables: Error allocating memory (u2d_destag)' + u2d_destag = 0._rk + call fm%send_data('u2d_destag', u2d_destag(_2D_FO_)) + end if + + if (v2d_destag_used) then + allocate(v2d_destag(E2DFIELD),stat=rc) + if (rc /= 0) stop 'finalize_register_processed_variables: Error allocating memory (v2d_destag)' + v2d_destag = 0._rk + call fm%send_data('v2d_destag', v2d_destag(_2D_FO_)) + end if + +#ifndef NO_BAROCLINIC + if (process_TEFmean) then + if (allocated(h_s_mean )) call fm%send_data('h_s_mean' , h_s_mean (_2D_FO_,:)) + if (allocated(hS_s_mean )) call fm%send_data('hS_s_mean' , hS_s_mean (_2D_FO_,:)) + if (allocated(hS2_s_mean )) call fm%send_data('hS2_s_mean' , hS2_s_mean (_2D_FO_,:)) + if (allocated(uu_s_mean )) call fm%send_data('uu_s_mean' , uu_s_mean (_2D_FO_,:)) + if (allocated(vv_s_mean )) call fm%send_data('vv_s_mean' , vv_s_mean (_2D_FO_,:)) + if (allocated(Sfluxu_s_mean)) call fm%send_data('Sfluxu_s_mean', Sfluxu_s_mean(_2D_FO_,:)) + if (allocated(Sfluxv_s_mean)) call fm%send_data('Sfluxv_s_mean', Sfluxv_s_mean(_2D_FO_,:)) + if (allocated(hpmS_s_mean )) call fm%send_data('hpmS_s_mean' , hpmS_s_mean (_2D_FO_,:)) + if (allocated(hnmS_s_mean )) call fm%send_data('hnmS_s_mean' , hnmS_s_mean (_2D_FO_,:)) + + if (allocated(fwf_s_int)) call fm%send_data('fwf_s_int', fwf_s_int (_2D_FO_,:)) + if (allocated(phymix_S_fwf_s_mean)) call fm%send_data('phymix_S_fwf_s_mean', phymix_S_fwf_s_mean(_2D_FO_,:)) + if (allocated(rvol_s_int)) call fm%send_data('rvol_s_int', rvol_s_int(_2D_FO_,:)) + if (allocated(phymix_S_riv_s_mean)) call fm%send_data('phymix_S_riv_s_mean', phymix_S_riv_s_mean(_2D_FO_,:)) + end if +#endif + return end subroutine finalize_register_processed_variables !EOC @@ -197,26 +401,52 @@ ! !IROUTINE: do_output_processing - read required variables ! ! !INTERFACE: - subroutine do_output_processing + subroutine do_output_processing(loop) ! ! !USES: use domain, only: az, au, av use variables_2d, only: z,D use variables_2d, only: U,V,DU,DV #ifndef NO_3D + use m3d, only: M use variables_3d, only: kmin,hn,uu,hun,vv,hvn +#ifndef NO_BAROCLINIC + use variables_3d, only: h_s, hS_s, hS2_s + use variables_3d, only: uu_s, vv_s + use variables_3d, only: Sfluxu_s, Sfluxv_s + use variables_3d, only: hpmS_s, hnmS_s + use variables_3d, only: fwf_s, phymix_S_fwf_s + use variables_3d, only: rvol_s, phymix_S_riv_s +#endif #endif IMPLICIT NONE ! ! !DESCRIPTION: ! ! !INPUT PARAMETERS: + integer, intent(in) :: loop ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding ! ! !LOCAL VARIABLES: REALTYPE, parameter :: vel_missing =-9999. + logical :: do_3d +#ifndef NO_BAROCLINIC + integer, save :: n_h_s=0 + integer, save :: n_hS_s=0 + integer, save :: n_hS2_s=0 + integer, save :: n_uu_s=0 + integer, save :: n_vv_s=0 + integer, save :: n_Sfluxu_s=0 + integer, save :: n_Sfluxv_s=0 + integer, save :: n_hpmS_s=0 + integer, save :: n_hnmS_s=0 + integer, save :: n_fwf_s=0 + integer, save :: n_phymix_S_fwf_s=0 + integer, save :: n_rvol_s=0 + integer, save :: n_phymix_S_riv_s=0 +#endif !EOP !------------------------------------------------------------------------- !BOC @@ -253,7 +483,7 @@ #endif - if (u2d_destag_use .and. v2d_destag_use) then + if (u2d_destag_now .and. v2d_destag_now) then call to_2d_u(imin,jmin,imax,jmax,az,U,DU,vel_missing, & imin,jmin,imax,jmax,u2d_destag) call to_2d_v(imin,jmin,imax,jmax,az,V,DV,vel_missing, & @@ -269,9 +499,57 @@ hvn,vv,vel_missing,v3d_destag) end if #endif +#endif + +#ifndef NO_3D + do_3d = (mod(loop,M) .eq. 0) +#ifndef NO_BAROCLINIC + if (process_TEFmean) then + call calc_mean(h_s , h_s_mean_now , n_h_s , h_s_mean ) + call calc_mean(hS_s , hS_s_mean_now , n_hS_s , hS_s_mean ) + call calc_mean(hS2_s , hS2_s_mean_now , n_hS2_s , hS2_s_mean ) + call calc_mean(uu_s , uu_s_mean_now , n_uu_s , uu_s_mean ) + call calc_mean(vv_s , vv_s_mean_now , n_vv_s , vv_s_mean ) + call calc_mean(Sfluxu_s, Sfluxu_s_mean_now, n_Sfluxu_s, Sfluxu_s_mean) + call calc_mean(Sfluxv_s, Sfluxv_s_mean_now, n_Sfluxv_s, Sfluxv_s_mean) + call calc_mean(hpmS_s , hpmS_s_mean_now , n_hpmS_s , hpmS_s_mean ) + call calc_mean(hnmS_s , hnmS_s_mean_now , n_hnmS_s , hnmS_s_mean ) + + call calc_mean(fwf_s, fwf_s_int_now, n_fwf_s , fwf_s_int, no_avg=.true.) + call calc_mean(phymix_S_fwf_s, phymix_S_fwf_s_mean_now, n_phymix_S_fwf_s, phymix_S_fwf_s_mean) + call calc_mean(rvol_s, rvol_s_int_now , n_rvol_s, rvol_s_int, no_avg=.true.) + call calc_mean(phymix_S_riv_s, phymix_S_riv_s_mean_now, n_phymix_S_riv_s, phymix_S_riv_s_mean) + end if +#endif #endif return + + contains +#ifndef NO_3D + subroutine calc_mean(snapshot, now, n, mean, no_avg) + REALTYPE, allocatable, intent(in) :: snapshot(:,:,:) + logical , intent(in) :: now + integer , intent(inout) :: n + REALTYPE, allocatable, intent(inout) :: mean(:,:,:) + logical , optional , intent(in) :: no_avg + logical :: do_avg + if (allocated(mean)) then + if (n .eq. 0) mean = _ZERO_ + if (do_3d) then + mean = mean + snapshot + n = n + 1 + end if + if (now) then + do_avg = .true. + if (present(no_avg)) do_avg = (.not. no_avg) + if (do_avg .and. n.gt.0) mean = mean / n + n = 0 + end if + end if + end subroutine calc_mean +#endif + end subroutine do_output_processing !EOC