From b98e72b234b083d6e3a75bb02228297eb521c59e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sat, 3 Sep 2016 11:59:35 -0400 Subject: [PATCH 1/4] Map tracers and vel. into depth space under an ice shelf This is a work around to map tracers and vel. into depth space under an ice shelf. The base of the ice shelf is determined based on the SSH field. The total water column is then dilated using the SSH field and the ocean depth. --- src/core/MOM.F90 | 2 +- src/diagnostics/MOM_diag_to_Z.F90 | 295 ++++++++++++++++++------------ 2 files changed, 177 insertions(+), 120 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2ae13af176..772c90619a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1159,7 +1159,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 From c626e645303db62d9f80be6ceb428e41d1857951 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Wed, 7 Sep 2016 10:39:14 -0400 Subject: [PATCH 2/4] Version of EPBL with iteration over OBL depth ready for testing. - Enable with flag: USE_MLD_ITERATION = .true. - Otherwise, old version of EPBL is employed. --- .../vertical/MOM_energetic_PBL.F90 | 281 ++++++++++++++++-- 1 file changed, 263 insertions(+), 18 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 95cdea0c70..09896fbe63 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,55 @@ 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. + LOGICAL :: Use_MLD_ITERATION=.false.!False to use old ePBL method. + 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. 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=50 ! Set maximum number of iterations. Probably + ! best as an input parameter, but will then need + ! to create allocatable arrays if storing + ! guess/found (as diagnostic); skipping for now. + real, dimension(SZK_(GV)+1) :: Vstar_Used, & ! 1D arrays used to store + Mixing_Length_Used ! Vstar and Mixing_Length + real, dimension(SZK_(GV)) :: TEMPdepth + !/BGR - remaining variables are related to tracking iteration statistics. + logical :: OBL_IT_STATS=.true. ! Flag for computing OBL iteration statistics + REAL :: ITguess(50), ITresult(50),ITmax(50),ITmin(50) ! Flag for storing guess/result + 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 +416,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 +443,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 +485,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 +502,40 @@ 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; + TEMPdepth(k)=max_MLD + 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 +582,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)**2 + 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 +791,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 +824,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 +861,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 +886,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 +1049,93 @@ 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 + ITmin(obl_it)=min_MLD + ITguess(obl_it) = MLD_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*,'Depth: ',TEMPdepth(nz) + !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 +1181,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 +1441,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 +1470,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-2') call get_param(param_file, mod, "ENABLE_THERMODYNAMICS", use_temperature, & "If true, temperature and salinity are used as state \n"//& @@ -1259,7 +1492,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 +1511,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 +1519,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) From 9df5b43b3a9c2a426c86ffdcb54d66d58b7190ef Mon Sep 17 00:00:00 2001 From: John Krasting Date: Thu, 8 Sep 2016 08:32:12 -0400 Subject: [PATCH 3/4] Mods for generic abiotic tracer routine - Compute sosga and pass to generic_tracer_source() since routine normalizes global alkalinity by the ratio of local to global SSS. - Add model time to generic_tracer_coupler_set call - Wrap MAX_FIELDS in ifndef block so that it can be modified via CPPDEFS - Include block from Jasmin John to reset the COBALT 'cased' tracer to 0. for all layers k>=2 --- config_src/dynamic/MOM_memory.h | 2 ++ src/tracer/MOM_generic_tracer.F90 | 29 ++++++++++++++++++++++++++--- 2 files changed, 28 insertions(+), 3 deletions(-) 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/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index e551d31bf5..f222a5ff59 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -60,6 +60,7 @@ module MOM_generic_tracer use MOM_grid, only : ocean_grid_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 @@ -373,6 +374,15 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, OBC 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 @@ -536,6 +546,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,G%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 @@ -605,11 +618,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) @@ -742,7 +760,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 @@ -846,6 +864,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 @@ -854,12 +874,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 MOM 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)//& From 0d934bb7346e400b0a2e34c053a32eac171c316e Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Thu, 8 Sep 2016 11:07:02 -0400 Subject: [PATCH 4/4] Updates to MOM_energetic_PBL iteration scheme. --- .../vertical/MOM_energetic_PBL.F90 | 43 ++++++++++--------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 09896fbe63..be1da1087b 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -345,7 +345,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !---------------------------------------------------------------------- !/BGR added Aug24,2016 for adding iteration to get boundary layer depth ! - needed to compute new mixing length. - LOGICAL :: Use_MLD_ITERATION=.false.!False to use old ePBL method. 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 @@ -356,7 +355,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! 3. Based on result adjusts the Max/Min ! and searches through the water column. ! - If using an accurate guess the iteration - ! is very quick. Otherwise it takes 5-10 + ! 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 @@ -368,21 +368,26 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! 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 + 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=50 ! Set maximum number of iterations. Probably - ! best as an input parameter, but will then need - ! to create allocatable arrays if storing + 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 - real, dimension(SZK_(GV)) :: TEMPdepth !/BGR - remaining variables are related to tracking iteration statistics. - logical :: OBL_IT_STATS=.true. ! Flag for computing OBL iteration statistics - REAL :: ITguess(50), ITresult(50),ITmax(50),ITmin(50) ! Flag for storing guess/result + 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) @@ -508,7 +513,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & 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; - TEMPdepth(k)=max_MLD enddo MIN_MLD = 0.0 !MIN_MLD will initialize as 0. @@ -531,8 +535,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & 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 + 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. !} @@ -590,7 +594,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & 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)**2 + 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 @@ -1051,9 +1055,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & 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 - ITmin(obl_it)=min_MLD - ITguess(obl_it) = MLD_GUESS + 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 @@ -1117,11 +1121,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !enddo !Activate to print out some output when not converged !{ - !print*,'Depth: ',TEMPdepth(nz) !print*,'Min/Max: ',ITmin(50),ITmax(50) !print*,'Guess/result: ',ITguess(50),ITresult(50) - print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,& - real(NOTCONVERGED)/real(CONVERGED) + !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,& + ! real(NOTCONVERGED)/real(CONVERGED) !} !stop !Kill if not converged during testing. else @@ -1473,7 +1476,7 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) 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-2') + 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"//&