This subroutine provides an interface for a user to use to modify the main code to alter the diffusivities as needed. The specific example implemented here augments the diffusivity for a specified range of latitude and coordinate potential density.
48 type(ocean_grid_type),
intent(in) :: g
49 type(verticalgrid_type),
intent(in) :: gv
50 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
51 type(thermo_var_ptrs),
intent(in) :: tv
54 type(unit_scale_type),
intent(in) :: us
55 type(user_change_diff_cs),
pointer :: cs
56 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(inout) :: kd_lay
58 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(inout) :: kd_int
60 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: t_f
62 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: s_f
64 real,
dimension(:,:,:),
optional,
pointer :: kd_int_add
68 real :: rcv(szi_(g),szk_(g))
69 real :: p_ref(szi_(g))
74 logical :: store_kd_add
75 integer,
dimension(2) :: eosdom
76 integer :: i, j, k, is, ie, js, je, nz
77 integer :: isd, ied, jsd, jed
81 character(len=200) :: mesg
83 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
84 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
86 if (.not.
associated(cs))
call mom_error(fatal,
"user_set_diffusivity: "//&
87 "Module must be initialized before it is used.")
89 use_eos =
associated(tv%eqn_of_state)
90 if (.not.use_eos)
return 91 store_kd_add = .false.
92 if (
present(kd_int_add)) store_kd_add =
associated(kd_int_add)
94 if (.not.range_ok(cs%lat_range))
then 95 write(mesg,
'(4(1pe15.6))') cs%lat_range(1:4)
96 call mom_error(fatal,
"user_set_diffusivity: bad latitude range: \n "//&
99 if (.not.range_ok(cs%rho_range))
then 100 write(mesg,
'(4(1pe15.6))') cs%rho_range(1:4)
101 call mom_error(fatal,
"user_set_diffusivity: bad density range: \n "//&
105 if (store_kd_add) kd_int_add(:,:,:) = 0.0
107 do i=is,ie ; p_ref(i) = tv%P_Ref ;
enddo 108 eosdom(:) = eos_domain(g%HI)
110 if (
present(t_f) .and.
present(s_f))
then 112 call calculate_density(t_f(:,j,k), s_f(:,j,k), p_ref, rcv(:,k), tv%eqn_of_state, eosdom)
116 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rcv(:,k), tv%eqn_of_state, eosdom)
120 if (
present(kd_lay))
then 121 do k=1,nz ;
do i=is,ie
122 if (cs%use_abs_lat)
then 123 lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
125 lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
127 rho_fn = val_weights(rcv(i,k), cs%rho_range)
128 if (rho_fn * lat_fn > 0.0) &
129 kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add * rho_fn * lat_fn
132 if (
present(kd_int))
then 133 do k=2,nz ;
do i=is,ie
134 if (cs%use_abs_lat)
then 135 lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
137 lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
139 rho_fn = val_weights( 0.5*(rcv(i,k-1) + rcv(i,k)), cs%rho_range)
140 if (rho_fn * lat_fn > 0.0)
then 141 kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add * rho_fn * lat_fn
142 if (store_kd_add) kd_int_add(i,j,k) = cs%Kd_add * rho_fn * lat_fn