diff --git a/config_src/dynamic/MOM_memory.h b/config_src/dynamic/MOM_memory.h index 43f952acd6..b2773188de 100644 --- a/config_src/dynamic/MOM_memory.h +++ b/config_src/dynamic/MOM_memory.h @@ -34,7 +34,9 @@ ! NJPROC_ is the number of processors in the ! y-direction. +#ifndef MAX_FIELDS_ #define MAX_FIELDS_ 50 +#endif ! The maximum permitted number (each) of ! restart variables, time derivatives, etc. ! This is mostly used for the size of pointer diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4baee857da..d0387f5011 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1162,7 +1162,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (Time_local + set_time(int(0.5*dt_therm)) > CS%Z_diag_time) then call enable_averaging(real(time_type_to_real(CS%Z_diag_interval)), & CS%Z_diag_time, CS%diag) - call calculate_Z_diag_fields(u, v, h, CS%dt_trans, & + call calculate_Z_diag_fields(u, v, h, ssh, CS%dt_trans, & G, GV, CS%diag_to_Z_CSp) CS%Z_diag_time = CS%Z_diag_time + CS%Z_diag_interval call disable_averaging(CS%diag) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 4c019bd9d8..9c0a1744bb 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -42,7 +42,7 @@ module MOM_diag_to_Z !* The boundaries always run through q grid points (x). * !* * !********+*********+*********+*********+*********+*********+*********+** - +use MOM_domains, only : pass_var use MOM_coms, only : reproducing_sum use MOM_diag_mediator, only : post_data, post_data_1d_k, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type, diag_axis_init @@ -162,27 +162,30 @@ function global_z_mean(var,G,CS,tracer) end function global_z_mean -subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) +subroutine calculate_Z_diag_fields(u, v, h, ssh_in, dt, G, GV, CS) type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_in real, intent(in) :: dt type(diag_to_Z_CS), pointer :: CS ! This subroutine maps tracers and velocities into depth space for diagnostics. ! Arguments: -! (in) u - zonal velocity component (m/s) -! (in) v - meridional velocity component (m/s) -! (in) h - layer thickness (meter or kg/m2) -! (in) dt - time difference in sec since last call to this subroutine -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - control structure returned by previous call to diagnostics_init +! (in) u - zonal velocity component (m/s) +! (in) v - meridional velocity component (m/s) +! (in) h - layer thickness (meter or kg/m2) +! (in) ssh_in - sea surface height (meter or kg/m2) +! (in) dt - time difference in sec since last call to this subroutine +! (in) G - ocean grid structure +! (in) GV - The ocean's vertical grid structure. +! (in) CS - control structure returned by previous call to diagnostics_init ! Note the deliberately reversed axes in h_f, u_f, v_f, and tr_f. + real :: ssh(SZI_(G),SZJ_(G)) ! copy of ssh_in (meter or kg/m2) real :: e(SZK_(G)+2) ! z-star interface heights (meter or kg/m2) real :: h_f(SZK_(G)+1,SZI_(G)) ! thicknesses of massive layers (meter or kg/m2) real :: u_f(SZK_(G)+1,SZIB_(G))! zonal velocity component in any massive layer @@ -191,7 +194,8 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) real :: tr_f(SZK_(G),max(CS%num_tr_used,1),SZI_(G)) ! tracer concentration in massive layers integer :: nk_valid(SZIB_(G)) ! number of massive layers in a column - real :: D_pt(SZIB_(G)) ! bottom depth (meter or kg/m2) + real :: D_pt(SZIB_(G)) ! bottom depth (meter or kg/m2) + real :: shelf_depth(SZIB_(G)) ! ice shelf depth (meter or kg/m2) real :: htot ! summed layer thicknesses (meter or kg/m2) real :: dilate ! proportion by which to dilate every layer real :: wt(SZK_(G)+1) ! fractional weight for each layer in the @@ -212,11 +216,15 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) integer :: k_top, k_bot, k_bot_prev integer :: i, j, k, k2, kz, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nk, m, nkml - + integer :: IsgB, IegB, JsgB, JegB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nk = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB nkml = max(GV%nkml, 1) Angstrom = GV%Angstrom + ssh(:,:) = ssh_in + ! Update the halos + call pass_var(ssh, G%Domain) if (.not.associated(CS)) call MOM_error(FATAL, & "diagnostic_fields_zstar: Module must be initialized before it is used.") @@ -231,10 +239,21 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) CS%u_z(I,j,kz) = CS%missing_vel enddo ; enddo ; enddo + do j=js,je ! Remove all massless layers. do I=Isq,Ieq - nk_valid(I) = 0 ; D_pt(I) = 0.5*(G%bathyT(i+1,j)+G%bathyT(i,j)) + nk_valid(I) = 0 + D_pt(I) = 0.5*(G%bathyT(i+1,j)+G%bathyT(i,j)) + ! GM: The following is workaround to make this subroutine works under an ice shelf + ! Zero ssh in the open ocean and get the depth of ice shelf (hence the abs below) + ! I am not sure if -0.1 is the best choice + shelf_depth(I) = 0.5*(ssh(i+1,j)+ssh(i,j)) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif enddo do k=1,nk ; do I=Isq,Ieq if ((G%mask2dCu(I,j) > 0.5) .and. (h(i,j,k)+h(i+1,j,k) > 4.0*Angstrom)) then @@ -247,13 +266,22 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! no-slip BBC in the output, if anything but piecewise constant is used. nk_valid(I) = nk_valid(I) + 1 ; k2 = nk_valid(I) h_f(k2,I) = Angstrom ; u_f(k2,I) = 0.0 + ! GM: D_pt is always slightly larger (by 1E-6 or so) than shelf_depth, so + ! I consider that the ice shelf is grounded when + ! shelf_depth(I) + 1.0E-3 > D_pt(i) + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 endif ; enddo + do I=Isq,Ieq ; if (nk_valid(I) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot*GV%H_to_m > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot*GV%H_to_m > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + !dilate = (D_pt(i) - 0.0)/htot + !write(*,*)MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -264,39 +292,44 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(I)) exit - if (linear_velocity_profiles) then - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(I)) .and. (k > nkml)) call & - find_limited_slope(u_f(:,I), e, slope, k) - endif - ! This is the piecewise linear form. - CS%u_z(I,j,kz) = wt(k) * (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) - enddo - if (k_bot > k_top) then ; k = k_bot - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(I)) .and. (k > nkml)) call & - find_limited_slope(u_f(:,I), e, slope, k) - ! This is the piecewise linear form. - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k) * & - (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - k_bot_prev = k_bot - else ! Use piecewise constant profiles. - CS%u_z(I,j,kz) = wt(k_top)*u_f(k_top,I) - do k=k_top+1,k_bot - CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) - enddo - endif ! linear profiles + !GM if top range that is being map is below the shelf, interpolate + ! otherwise keep missing_vel + if (CS%Z_int(kz)<=-shelf_depth(I)) then + + if (linear_velocity_profiles) then + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(I)) .and. (k > nkml)) call & + find_limited_slope(u_f(:,I), e, slope, k) + endif + ! This is the piecewise linear form. + CS%u_z(I,j,kz) = wt(k) * (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) + enddo + if (k_bot > k_top) then ; k = k_bot + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(I)) .and. (k > nkml)) call & + find_limited_slope(u_f(:,I), e, slope, k) + ! This is the piecewise linear form. + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k) * & + (u_f(k,I) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + k_bot_prev = k_bot + else ! Use piecewise constant profiles. + CS%u_z(I,j,kz) = wt(k_top)*u_f(k_top,I) + do k=k_top+1,k_bot + CS%u_z(I,j,kz) = CS%u_z(I,j,kz) + wt(k)*u_f(k,I) + enddo + endif ! linear profiles + endif ! below shelf enddo ! kz-loop endif ; enddo ! I-loop and mask enddo ! j-loop @@ -314,6 +347,12 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! Remove all massless layers. do i=is,ie nk_valid(i) = 0 ; D_pt(i) = 0.5*(G%bathyT(i,j)+G%bathyT(i,j+1)) + shelf_depth(I) = 0.5*(ssh(i,j)+ssh(i,j+1)) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif enddo do k=1,nk ; do i=is,ie if ((G%mask2dCv(i,j) > 0.5) .and. (h(i,j,k)+h(i,j+1,k) > 4.0*Angstrom)) then @@ -326,13 +365,16 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) ! no-slip BBC in the output, if anything but piecewise constant is used. nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i) h_f(k2,i) = Angstrom ; v_f(k2,i) = 0.0 + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 endif ; enddo do i=is,ie ; if (nk_valid(i) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -342,41 +384,44 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) call find_overlap(e, CS%Z_int(kz), CS%Z_int(kz+1), nk_valid(i), & k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(i)) exit - - if (linear_velocity_profiles) then - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(v_f(:,i), e, slope, k) - endif - ! This is the piecewise linear form. - CS%v_z(i,J,kz) = wt(k) * (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) - enddo - if (k_bot > k_top) then ; k = k_bot - ! Calculate the intra-cell profile. - slope = 0.0 ! ; curv = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(v_f(:,i), e, slope, k) - ! This is the piecewise linear form. - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k) * & - (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - k_bot_prev = k_bot - else ! Use piecewise constant profiles. - CS%v_z(i,J,kz) = wt(k_top)*v_f(k_top,i) - do k=k_top+1,k_bot - CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) - enddo - endif ! linear profiles - enddo ! kz-loop + !GM if top range that is being map is below the shelf, interpolate + ! otherwise keep missing_vel + if (CS%Z_int(kz)<=-shelf_depth(I)) then + if (linear_velocity_profiles) then + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(v_f(:,i), e, slope, k) + endif + ! This is the piecewise linear form. + CS%v_z(i,J,kz) = wt(k) * (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) + enddo + if (k_bot > k_top) then ; k = k_bot + ! Calculate the intra-cell profile. + slope = 0.0 ! ; curv = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(v_f(:,i), e, slope, k) + ! This is the piecewise linear form. + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k) * & + (v_f(k,i) + 0.5*slope*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + k_bot_prev = k_bot + else ! Use piecewise constant profiles. + CS%v_z(i,J,kz) = wt(k_top)*v_f(k_top,i) + do k=k_top+1,k_bot + CS%v_z(i,J,kz) = CS%v_z(i,J,kz) + wt(k)*v_f(k,i) + enddo + endif ! linear profiles + endif ! below shelf + enddo ! kz-loop endif ; enddo ! i-loop and mask enddo ! J-loop @@ -392,11 +437,20 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) do j=js,je ! Remove all massless layers. - do i=is,ie ; nk_valid(i) = 0 ; D_pt(i) = G%bathyT(i,j) ; enddo + do i=is,ie + nk_valid(i) = 0 ; D_pt(i) = G%bathyT(i,j) + shelf_depth(I) = ssh(i,j) + if (shelf_depth(I) > -0.1) then ! open ocean + shelf_depth(I) = 0.0 + else ! ice shelf + shelf_depth(I) = abs(shelf_depth(I)) + endif + enddo do k=1,nk ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 2.0*Angstrom)) then nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i) h_f(k2,i) = h(i,j,k) + if (shelf_depth(I) + 1.0E-3 > D_pt(i)) nk_valid(I)=0 do m=1,CS%num_tr_used ; tr_f(k2,m,i) = CS%tr_model(m)%p(i,j,k) ; enddo endif enddo ; enddo @@ -404,8 +458,10 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) do i=is,ie ; if (nk_valid(i) > 0) then ! Calculate the z* interface heights for tracers. htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo - dilate = 0.0 ; if (htot > 2.0*Angstrom) dilate = (D_pt(i) - 0.0) / htot - + dilate = 0.0 + if (htot > 2.0*Angstrom) then + dilate = MAX((D_pt(i) - shelf_depth(i)),Angstrom)/htot + endif e(nk_valid(i)+1) = -D_pt(i) do k=nk_valid(i),1,-1 ; e(K) = e(K+1) + h_f(k,i)*dilate ; enddo @@ -415,38 +471,39 @@ subroutine calculate_Z_diag_fields(u, v, h, dt, G, GV, CS) call find_overlap(e, CS%Z_int(kz), CS%Z_int(kz+1), nk_valid(i), & k_bot, k_top, k_bot, wt, z1, z2) if (k_top>nk_valid(i)) exit - - do m=1,CS%num_tr_used - k = k_top - if (k /= k_bot_prev) then - ! Calculate the intra-cell profile. - sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) - endif - ! This is the piecewise linear form. - CS%tr_z(m)%p(i,j,kz) = wt(k) * & - (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - do k=k_top+1,k_bot-1 - CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k)*tr_f(k,m,i) - enddo - if (k_bot > k_top) then - k = k_bot - ! Calculate the intra-cell profile. - sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 - if ((k < nk_valid(i)) .and. (k > nkml)) call & - find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) - ! This is the piecewise linear form. - CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k) * & - (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) - ! For the piecewise parabolic form add the following... - ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) - endif - enddo - k_bot_prev = k_bot - enddo ! kz-loop + if (CS%Z_int(kz)<=-shelf_depth(I)) then + do m=1,CS%num_tr_used + k = k_top + if (k /= k_bot_prev) then + ! Calculate the intra-cell profile. + sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) + endif + ! This is the piecewise linear form. + CS%tr_z(m)%p(i,j,kz) = wt(k) * & + (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + do k=k_top+1,k_bot-1 + CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k)*tr_f(k,m,i) + enddo + if (k_bot > k_top) then + k = k_bot + ! Calculate the intra-cell profile. + sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0 + if ((k < nk_valid(i)) .and. (k > nkml)) call & + find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k) + ! This is the piecewise linear form. + CS%tr_z(m)%p(i,j,kz) = CS%tr_z(m)%p(i,j,kz) + wt(k) * & + (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k))) + ! For the piecewise parabolic form add the following... + ! + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2)) + endif + enddo + k_bot_prev = k_bot + endif ! below shelf + enddo ! kz-loop endif ; enddo ! i-loop and mask enddo ! j-loop diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 95cdea0c70..be1da1087b 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -116,13 +116,16 @@ module MOM_energetic_PBL ! diffusivity in the planetary boundary layer. type(time_type), pointer :: Time ! A pointer to the ocean model's clock. logical :: TKE_diagnostics = .false. - + LOGICAL :: Use_MLD_ITERATION=.false.!False to use old ePBL method. + logical :: Mixing_Diagnostics = .false. ! Will be true when outputing mixing + ! length and velocity scale type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the ! timing of diagnostic output. ! These are terms in the mixed layer TKE budget, all in J m-2 = kg s-2. real, allocatable, dimension(:,:) :: & ML_depth, & ! The mixed layer depth in m. + ML_depth2, & ! The mixed layer depth in m. diag_TKE_wind, & ! The wind source of TKE. diag_TKE_MKE, & ! The resolved KE source of TKE. diag_TKE_conv, & ! The convective source of TKE. @@ -132,10 +135,14 @@ module MOM_energetic_PBL diag_TKE_conv_decay, & ! The decay of convective TKE. diag_TKE_mixing ! The work done by TKE to deepen ! the mixed layer. + real, allocatable, dimension(:,:,:) :: & + Velocity_Scale, & ! The velocity scale used in getting Kd + Mixing_Length ! The length scale used in getting Kd integer :: id_ML_depth = -1, id_TKE_wind = -1, id_TKE_mixing = -1 integer :: id_TKE_MKE = -1, id_TKE_conv = -1, id_TKE_forcing = -1 integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1 integer :: id_Hsfc_used = -1 + integer :: id_Mixing_Length = -1, id_Velocity_Scale = -1 end type energetic_PBL_CS integer :: num_msg = 0, max_msg = 2 @@ -335,7 +342,60 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & logical :: write_diags ! If true, write out diagnostics with this step. logical :: reset_diags ! If true, zero out the accumulated diagnostics. ! detrainment, in units of m. - + !---------------------------------------------------------------------- + !/BGR added Aug24,2016 for adding iteration to get boundary layer depth + ! - needed to compute new mixing length. + real :: MLD_GUESS, MLD_FOUND ! Mixing Layer depth guessed/found for iteration + real :: MAX_MLD, MIN_MLD ! Iteration bounds which are adjusted at each step + ! - These are initialized based on surface/bottom + ! 1. The iteration guesses a value (possibly from + ! prev step or neighbor). + ! 2. The iteration checks if value is converged, + ! too shallow, or too deep. + ! 3. Based on result adjusts the Max/Min + ! and searches through the water column. + ! - If using an accurate guess the iteration + ! is very quick (e.g. if MLD doesn't change + ! over timestep). Otherwise it takes 5-10 + ! passes, but has a high convergence rate. + ! Other iteration may be tried, but this + ! method seems to rarely fail and the added + ! cost is likely not significant. Additionally, + ! when it fails it does so in a reasonable + ! manner giving a usable guess. When it + ! does fail, it is due to convection within + ! the boundary. Likely, a new method e.g. + ! surface_disconnect, can improve this. + logical :: FIRST_OBL ! Flag for computing "found" Mixing layer depth + logical :: OBL_CONVERGED ! Flag for convergence of MLD + INTEGER :: OBL_IT ! Iteration counter + REAL :: MinMixLen=1.0 ! Minimum value for mixing length (must be non-zero + ! for iteration to converge when OBL shoals). + ! At present, this is the value used for convection + ! in the ocean interior. + INTEGER :: MAX_OBL_IT=20 ! Set maximum number of iterations. Probably + ! best as an input parameter, but then may want + ! to use allocatable arrays if storing + ! guess/found (as diagnostic); skipping for now. + ! In reality, the maximum number of guesses + ! needed is set by: + ! DEPTH/2^M < DZ + ! where M is the number of guesses + ! e.g. M=12 for DEPTH=4000m and DZ=1m + real, dimension(SZK_(GV)+1) :: Vstar_Used, & ! 1D arrays used to store + Mixing_Length_Used ! Vstar and Mixing_Length + !/BGR - remaining variables are related to tracking iteration statistics. + logical :: OBL_IT_STATS=.false. ! Flag for computing OBL iteration statistics + REAL :: ITguess(20), ITresult(20),ITmax(20),ITmin(20) ! Flag for storing guess/result + ! should have dim=MAX_OBL_IT + integer, save :: MAXIT=0 ! Stores maximum number of iterations + integer, save :: MINIT=1e8 ! Stokes minimum number of iterations + integer, save :: SUMIT=0 ! Stores total iterations (summed over all) + integer, save :: NUMIT=0 ! Stores number of times iterated + !e.g. Average iterations = SUMIT/NUMIT + integer, save :: CONVERGED! + integer, save :: NOTCONVERGED! + !-End BGR iteration parameters----------------------------------------- logical :: debug=.false. ! Change this hard-coded value for debugging. ! The following arrays are used only for debugging purposes. real :: dPE_debug, mixing_debug @@ -361,6 +421,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & h_neglect = GV%H_subroundoff + if(.not.CS%Use_MLD_Iteration) MAX_OBL_IT=1 C1_3 = 1.0 / 3.0 dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag IdtdR0 = 1.0 / (dt__diag * GV%Rho0) @@ -387,6 +448,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & CS%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced_forcing(i,j) = 0.0 enddo ; enddo endif + IF (CS%Mixing_Diagnostics) then + CS%Mixing_Length(:,:,:)=0.0 + CS%Velocity_Scale(:,:,:)=0.0 + endif !!OMP end parallel endif @@ -425,7 +490,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & T(i,k) = tv%T(i,j,k) ; S(i,k) = tv%S(i,j,k) Kd(i,K) = 0.0 enddo ; enddo - do i=is,ie ; CS%ML_depth(i,j) = h(i,1)*GV%H_to_m ; sfc_connected(i) = .true. ; enddo + do i=is,ie + CS%ML_depth(i,j) = h(i,1)*GV%H_to_m ; + !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m ; + sfc_connected(i) = .true. ; + enddo if (debug) then mech_TKE_k(:,:) = 0.0 ; conv_PErel_k(:,:) = 0.0 @@ -438,6 +507,39 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! and ustar and wstar available to drive mixing at the first interior ! interface. do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + !/The following lines are for the iteration over MLD + !{ + OBL_CONVERGED=.false.!Initialize convergence state to 'false' + MAX_MLD = 0.0;!MAX_MLD will initialized as ocean bottom depth + do k=1,nz ; + MAX_MLD = MAX_MLD + h(i,k)*GV%H_to_m; + enddo + MIN_MLD = 0.0 !MIN_MLD will initialize as 0. + + !/BGR: Add MLD_GUESS based on stored previous value. + ! note that this is different from ML_Depth already + ! computed by EPBL, need to figure out why. + if (CS%ML_Depth2(i,j).gt.1.) then + !If prev value is present use for guess. + MLD_GUESS=CS%ML_Depth2(i,j) + else + !Otherwise guess middle of water column. + MLD_GUESS = 0.5 * (MIN_MLD+MAX_MLD) + endif + !/BGR: May add user-input bounds for max/min MLD + + DO OBL_IT=1,MAX_OBL_IT ! Iterate up to MAX_OBL_IT times + IF (.not. OBL_CONVERGED) THEN !Cycle through EPBL if not converged + CS%ML_depth(i,j) = h(i,1)*GV%H_to_m ; + !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m ; + sfc_connected(i) = .true. ; + ! Store in 1D arrays cleared out each iteration. Only write in + ! 3D arrays after convergence. + VSTAR_USED(:)=0.0 + MIXING_LENGTH_USED(:)=0.0 + IF (.not.CS%Use_MLD_Iteration) OBL_CONVERGED=.true. + !/This ends the new code for the iteration. + !} U_Star = fluxes%ustar(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & @@ -484,15 +586,28 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & nstar_k(:) = 0.0 ; nstar_k(1) = CS%nstar endif - h_sum(i) = H_neglect - do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo - I_hs = 0.0 ; if (h_sum(i) > 0.0) I_hs = 1.0 / h_sum(i) + if (.not.CS%Use_MLD_Iteration) then + h_sum(i) = H_neglect + do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo + I_hs = 0.0 ; if (h_sum(i) > 0.0) I_hs = 1.0 / h_sum(i) - h_bot(i) = 0.0 ; hb_hs(i,nz+1) = 0.0 - do k=nz,1,-1 - h_bot(i) = h_bot(i) + h(i,k) - hb_hs(i,K) = GV%H_to_m * (h_bot(i)*I_hs) - enddo + h_bot(i) = 0.0 ; hb_hs(i,nz+1) = 0.0 + do k=nz,1,-1 + h_bot(i) = h_bot(i) + h(i,k) + hb_hs(i,K) = GV%H_to_m * (h_bot(i)*I_hs) + enddo + else + !New method where mixing length is reduced based on MLD + I_hs=1.0/MLD_GUESS + h_sum(i) = 0.0 + h_bot(i) = 0.0 ; hb_hs(i,1:nz+1) = 0.0 + do k=2,nz + h_sum(i)=h_sum(i)+h(i,k) + h_bot(i) = max(0.0,MLD_GUESS - h_sum(i)) + hb_hs(i,K) = GV%H_to_m * (h_bot(i)*I_hs)**2!Notice the square + ! makes KPP-like. + enddo + endif ! endif ; enddo ! Note the outer i-loop and inner k-loop loop order!!! @@ -680,11 +795,20 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*conv_PErel(i) if (TKE_here > 0.0) then vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 - Kd_guess0 = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & + Mixing_Length_Used(k) = MAX(MinMixLen,((h_tt*hb_hs(i,K))*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar)) + !Note setting Kd_guess0 to Mixing_Length_Used(K) here will + ! change the answers. Therefore, skipping that. + if (.not.CS%Use_MLD_Iteration) then + Kd_guess0 = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar) + else + Kd_guess0 = vstar * vonKar * Mixing_Length_Used(k) + endif else vstar = 0.0 ; Kd_guess0 = 0.0 endif + Vstar_Used(k) = vstar!Track vstar Kddt_h_g0 = Kd_guess0*dt_h call find_PE_chg(Kddt_h_g0, h(i,k), b_den_1(i), dTe_term, dSe_term, & @@ -704,11 +828,20 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & TKE_here = mech_TKE(i) + CS%wstar_ustar_coef*(conv_PErel(i)-PE_chg_max) if (TKE_here > 0.0) then vstar = CS%vstar_scale_fac * (I_dtrho*TKE_here)**C1_3 - Kd(i,k) = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & - ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar) + Mixing_Length_Used(k) = max(MinMixLen,((h_tt*hb_hs(i,K))*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar)) + if (.not.CS%Use_MLD_Iteration) then + ! Note again (as prev) that using Mixing_Length_Used here + ! instead of redoing the computation will change answers... + Kd(i,k) = vstar * vonKar * ((h_tt*hb_hs(i,K))*vstar) / & + ((CS%Ekman_scale_coef * absf(i)) * (h_tt*hb_hs(i,K)) + vstar) + else + Kd(i,k) = vstar * vonKar * Mixing_Length_Used(k) + endif else vstar = 0.0 ; Kd(i,k) = 0.0 endif + Vstar_Used(k) = vstar call find_PE_chg(Kd(i,k)*dt_h, h(i,k), b_den_1(i), dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE(i,k), dS_to_dPE(i,k), & @@ -732,8 +865,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & CS%diag_TKE_conv(i,j) = CS%diag_TKE_conv(i,j) - CS%nstar*dPE_conv * IdtdR0 CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + MKE_src * IdtdR0 endif - if (sfc_connected(i)) CS%ML_depth(i,J) = CS%ML_depth(i,J) + & - GV%H_to_m * h(i,k) + if (sfc_connected(i)) then + CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) + !CS%ML_depth2(i,j) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) + endif Kddt_h(K) = Kd(i,k)*dt_h elseif (tot_TKE + (MKE_src - PE_chg_g0) >= 0.0) then @@ -755,8 +890,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & tot_TKE = TKE_reduc*tot_TKE mech_TKE(i) = TKE_reduc*(mech_TKE(i) + MKE_src) conv_PErel(i) = TKE_reduc*conv_PErel(i) - if (sfc_connected(i)) CS%ML_depth(i,J) = CS%ML_depth(i,J) + & - GV%H_to_m * h(i,k) + if (sfc_connected(i)) then + CS%ML_depth(i,J) = CS%ML_depth(i,J) + GV%H_to_m * h(i,k) + !CS%ML_depth2(i,J) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k) + endif elseif (tot_TKE == 0.0) then ! This can arise if nstar_FC = 0. Kd(i,k) = 0.0 ; Kddt_h(K) = 0.0 @@ -916,6 +1053,92 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & mixing_debug = dPE_debug * IdtdR0 endif k = nz ! This is here to allow a breakpoint to be set. + + !/BGR: The following lines are used for the iteration + ITmax(obl_it)=max_MLD ! Track max } + ITmin(obl_it)=min_MLD ! Track min } For debug purpose + ITguess(obl_it) = MLD_GUESS ! Track guess } + MLD_FOUND=0.0;FIRST_OBL=.true.; + ! MLD_FOUND=CS%ML_depth2(i,J) + do k=2,nz + IF (FIRST_OBL) then !Breaks when OBL found + if (Vstar_Used(k).gt.1.e-10 .and. k.lt.nz) then + !If within OBL, keep integrating to find OBL + !/ Add thickness of level above to OBL depth. + MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m + !1. Check if guess was too shallow + !Adding -1 m as cushion, will help avoid + ! non-convergence flag when nearly converged. + elseif (MLD_FOUND-1.0.GT.MLD_GUESS) then + !elseif (MLD_FOUND.GT.MLD_GUESS) then + !/ Guess was too shallow, set new minimum guess + MIN_MLD=MLD_GUESS + FIRST_OBL=.false. !Break OBL loop + !2. Check if Guess minus found MLD + ! is less than thickness of level (=converged) + ! - We could try to add a more precise + ! value for found MLD, but seems difficult to + ! to contrain beyond within a level. + elseif ((MLD_GUESS-MLD_FOUND).lt.max(1.,h(i,k-1)*GV%H_to_m)) then + ! Converged. Exit iteration. + FIRST_OBL=.false.!Break OBL loop + OBL_CONVERGED=.true.!Break convergence loop + ! Testing Output, comment for use. + !print*,'Converged--------' + !print*,MLD_FOUND,MLD_GUESS + !/ + IF (OBL_IT_STATS) then !Compute iteration statistics + MAXIT=max(MAXIT,obl_it) + MINIT=min(MINIT,obl_it) + SUMIT=SUMIT+obl_it + NUMIT=NUMIT+1 + print*,MAXIT,MINIT,SUMIT/NUMIT + ENDIF + !BGR this can be where MLD is stored for next time... + CS%ML_Depth2(i,j)=MLD_GUESS + !/ + !2. If not, guess was too deep + else + !Guess was too deep, set new maximum guess + MAX_MLD=MLD_GUESS !We know this guess was too deep + FIRST_OBL=.false.!Break OBL loop + endif + endif + enddo + ! For next pass, guess average of minimum and maximum values. + MLD_GUESS = MIN_MLD*0.5 + MAX_MLD*0.5 + ITresult(obl_it)=MLD_FOUND + ENDIF!Check for convergence loop + ENDDO!Iteration loop + if (.not.OBL_CONVERGED) then + !/Temp output, warn that EPBL didn't converge + !/Print guess/found for every iteration step + !print*,'EPBL MLD DID NOT CONVERGE' + NOTCONVERGED=NOTCONVERGED+1 + !do obl_it=1,max_obl_it + ! print*,ITmin(obl_it),ITmax(obl_it) + ! print*,obl_it,ITguess(obl_it),ITresult(obl_it) + !enddo + !Activate to print out some output when not converged + !{ + !print*,'Min/Max: ',ITmin(50),ITmax(50) + !print*,'Guess/result: ',ITguess(50),ITresult(50) + !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,& + ! real(NOTCONVERGED)/real(CONVERGED) + !} + !stop !Kill if not converged during testing. + else + CONVERGED=CONVERGED+1 + endif + if (cs%Mixing_Diagnostics) then + !Write to 3-D for outputing Mixing length and + ! velocity scale. + do k=1,nz + cs%Mixing_Length(i,j,k)=Mixing_Length_Used(k) + cs%Velocity_Scale(i,j,k)=Vstar_Used(k) + enddo + endif + else ! For masked points, Kd_int must still be set (to 0) because it has intent(out). do K=1,nz+1 @@ -961,6 +1184,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & call post_data(CS%id_TKE_conv_decay, CS%diag_TKE_conv_decay, CS%diag) if (CS%id_Hsfc_used > 0) & call post_data(CS%id_Hsfc_used, Hsfc_used, CS%diag) + if (CS%id_Mixing_Length > 0) & + call post_data(CS%id_Mixing_Length, CS%Mixing_Length, CS%diag) + if (CS%id_Velocity_Scale >0) & + call post_data(CS%id_Velocity_Scale, CS%Velocity_Scale, CS%diag) endif end subroutine energetic_PBL @@ -1217,6 +1444,11 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) "of the diffusive length scale by rotation. Making this larger\n"//& "decreases the PBL diffusivity.", & "units=nondim", default=1.0) + call get_param(param_file, mod, "USE_MLD_ITERATION", CS%USE_MLD_ITERATION, & + "A logical that determines whether or not to use the\n"//& + "MLD to set the EPBL length scale.", default=.false.) + + ! This gives a minimum decay scale that is typically much less than Angstrom. CS%ustar_min = 2e-4*CS%omega*(GV%Angstrom_z + GV%H_to_m*GV%H_subroundoff) @@ -1241,6 +1473,10 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) Time, 'Convective energy decay sink of mixed layer TKE', 'meter3 second-3') CS%id_Hsfc_used = register_diag_field('ocean_model', 'ePBL_Hs_used', diag%axesT1, & Time, 'Surface region thickness that is used', 'meter') + CS%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, & + Time, 'Mixing Length that is used', 'meter') + CS%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, & + Time, 'Velocity Scale that is used.', 'meter second') call get_param(param_file, mod, "ENABLE_THERMODYNAMICS", use_temperature, & "If true, temperature and salinity are used as state \n"//& @@ -1259,7 +1495,16 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) CS%TKE_diagnostics = .true. endif + if (max(CS%id_Mixing_Length,CS%id_Velocity_Scale)>0) then + call safe_alloc_alloc(CS%Velocity_Scale,isd,ied,jsd,jed,SZK_(GV)+1) + call safe_alloc_alloc(CS%Mixing_Length,isd,ied,jsd,jed,SZK_(GV)+1) + CS%Velocity_Scale(:,:,:)=0.0 + CS%Mixing_Length(:,:,:)=0.0 + CS%mixing_diagnostics=.true. + endif call safe_alloc_alloc(CS%ML_depth, isd, ied, jsd, jed) + call safe_alloc_alloc(CS%ML_depth2, isd, ied, jsd, jed) + end subroutine energetic_PBL_init @@ -1269,6 +1514,7 @@ subroutine energetic_PBL_end(CS) if (.not.associated(CS)) return if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) + if (allocated(CS%ML_depth2)) deallocate(CS%ML_depth2) if (allocated(CS%diag_TKE_wind)) deallocate(CS%diag_TKE_wind) if (allocated(CS%diag_TKE_MKE)) deallocate(CS%diag_TKE_MKE) if (allocated(CS%diag_TKE_conv)) deallocate(CS%diag_TKE_conv) @@ -1276,6 +1522,8 @@ subroutine energetic_PBL_end(CS) if (allocated(CS%diag_TKE_mixing)) deallocate(CS%diag_TKE_mixing) if (allocated(CS%diag_TKE_mech_decay)) deallocate(CS%diag_TKE_mech_decay) if (allocated(CS%diag_TKE_conv_decay)) deallocate(CS%diag_TKE_conv_decay) + if (allocated(CS%Mixing_Length)) deallocate(CS%Mixing_Length) + if (allocated(CS%Velocity_Scale)) deallocate(CS%Velocity_Scale) deallocate(CS) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 6c2c6c2bef..d45cd78ca5 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -61,6 +61,7 @@ module MOM_generic_tracer use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS + use MOM_spatial_means, only : global_area_mean use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS use MOM_time_manager, only : time_type, get_time, set_time @@ -375,6 +376,15 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia endif enddo; enddo ; enddo + !jgj: Reset CASED to 0 below K=1 + if(trim(g_tracer_name) .eq. 'cased') then + do k=2,nk ; do j=jsc,jec ; do i=isc,iec + if(tr_ptr(i,j,k) .ne. CS%tracer_land_val) then + tr_ptr(i,j,k) = 0.0 + endif + enddo; enddo ; enddo + endif + else !Do it old way if the tracer is not registered to start from a specific source file. !This path should be deprecated if all generic tracers are required to start from specified sources. if (len_trim(CS%IC_file) > 0) then @@ -538,6 +548,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G character(len=fm_string_len) :: g_tracer_name real, dimension(:,:), pointer :: stf_array,trunoff_array,runoff_tracer_flux_array + real :: surface_field(SZI_(G),SZJ_(G)) + real :: sosga + real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke) :: rho_dzt, dzt real, dimension(G%isd:G%ied,G%jsd:G%jed) :: hblt_depth integer :: i, j, k, isc, iec, jsc, jec, nk @@ -607,11 +620,16 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G enddo; enddo ; enddo + do j=jsc,jec ; do i=isc,iec + surface_field(i,j) = tv%S(i,j,1) + enddo ; enddo + sosga = global_area_mean(surface_field, G) + ! !Calculate tendencies (i.e., field changes at dt) from the sources / sinks ! - call generic_tracer_source(tv%T,tv%S,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,& + call generic_tracer_source(tv%T,tv%S,sosga,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,& G%areaT,get_diag_time_end(CS%diag),& optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band) @@ -744,7 +762,7 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, yg real, dimension(:), intent(out) :: gmin,gmax real, dimension(:), intent(out) :: xgmin, ygmin, zgmin, xgmax, ygmax, zgmax type(ocean_grid_type), intent(in) :: G - type(MOM_generic_tracer_CS), pointer :: CS + type(MOM_generic_tracer_CS), pointer :: CS character(len=*), dimension(:), intent(out) :: names character(len=*), dimension(:), intent(out) :: units integer :: MOM_generic_tracer_min_max @@ -848,6 +866,8 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS) ! (in) CS - The control structure returned by a previous call to ! register_MOM_generic_tracer. + real :: sosga + character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_surface_state' real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke,1) :: rho0 type(g_tracer_type), pointer :: g_tracer @@ -856,12 +876,15 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS) !nnz: fake rho0 rho0=1.0 + sosga = global_area_mean(state%SSS, G) + call generic_tracer_coupler_set(state%tr_fields,& ST=state%SST,& SS=state%SSS,& + sosga=sosga, & rho=rho0,& !nnz: required for MOM5 and previous versions. ilb=G%isd, jlb=G%jsd,& - tau=1) + tau=1,model_time=get_diag_time_end(CS%diag)) !Output diagnostics via diag_manager for all tracers in this module ! if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&