This function initializes a tracer by reading a Z-space file, returning .true. if this appears to have been successful, and false otherwise.
31 logical :: tracer_Z_init
32 type(ocean_grid_type),
intent(in) :: G
33 type(unit_scale_type),
intent(in) :: US
34 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
36 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
38 character(len=*),
intent(in) :: filename
39 character(len=*),
intent(in) :: tr_name
41 real,
optional,
intent(in) :: missing_val
42 real,
optional,
intent(in) :: land_val
47 integer,
save :: init_calls = 0
49 #include "version_variable.h"
50 character(len=40) :: mdl =
"MOM_tracer_Z_init"
51 character(len=256) :: mesg
53 real,
allocatable,
dimension(:,:,:) :: &
55 real,
allocatable,
dimension(:) :: &
74 logical :: has_edges, use_missing, zero_surface
75 character(len=80) :: loc_msg
76 integer :: k_top, k_bot, k_bot_prev
77 integer :: i, j, k, kz, is, ie, js, je, nz, nz_in
78 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
80 landval = 0.0 ;
if (
present(land_val)) landval = land_val
82 zero_surface = .false.
85 if (
present(missing_val))
then
86 use_missing = .true. ; missing = missing_val
91 call read_z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, &
92 missing, scale=us%m_to_Z)
94 tracer_z_init = .false.
98 allocate(tr_in(g%isd:g%ied,g%jsd:g%jed,nz_in)) ; tr_in(:,:,:) = 0.0
99 allocate(tr_1d(nz_in)) ; tr_1d(:) = 0.0
100 call mom_read_data(filename, tr_name, tr_in(:,:,:), g%Domain)
104 if (
present(missing_val))
then
105 do j=js,je ;
do i=is,ie
106 if (g%mask2dT(i,j) == 0.0)
then
107 tr_in(i,j,1) = landval
108 elseif (abs(tr_in(i,j,1) - missing_val) <= 1e-6*abs(missing_val))
then
109 write(loc_msg,
'(f7.2," N ",f7.2," E")') g%geoLatT(i,j), g%geoLonT(i,j)
110 if (zero_surface)
then
111 call mom_error(warning,
"tracer_Z_init: Missing value of "// &
112 trim(tr_name)//
" found in an ocean point at "//trim(loc_msg)// &
113 " in "//trim(filename) )
116 call mom_error(fatal,
"tracer_Z_init: Missing value of "// &
117 trim(tr_name)//
" found in an ocean point at "//trim(loc_msg)// &
118 " in "//trim(filename) )
122 do k=2,nz_in ;
do j=js,je ;
do i=is,ie
123 if (abs(tr_in(i,j,k) - missing_val) <= 1e-6*abs(missing_val)) &
124 tr_in(i,j,k) = tr_in(i,j,k-1)
125 enddo ;
enddo ;
enddo
128 allocate(wt(nz_in+1)) ;
allocate(z1(nz_in+1)) ;
allocate(z2(nz_in+1))
134 do i=is,ie ; htot(i) = 0.0 ;
enddo
135 do k=1,nz ;
do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ;
enddo ;
enddo
137 do i=is,ie ;
if (g%mask2dT(i,j)*htot(i) > 0.0)
then
139 dilate = (g%bathyT(i,j) - 0.0) / htot(i)
140 e(nz+1) = -g%bathyT(i,j)
141 do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ;
enddo
144 do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ;
enddo
145 k_bot = 1 ; k_bot_prev = -1
147 if (e(k+1) > z_edges(1))
then
149 elseif (e(k) < z_edges(nz_in+1))
then
150 tr(i,j,k) = tr_1d(nz_in)
152 call find_overlap(z_edges, e(k), e(k+1), nz_in, &
153 k_bot, k_top, k_bot, wt, z1, z2)
155 if (kz /= k_bot_prev)
then
158 if ((kz < nz_in) .and. (kz > 1))
call &
159 find_limited_slope(tr_1d, z_edges, sl_tr, kz)
162 tr(i,j,k) = wt(kz) * &
163 (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
166 do kz=k_top+1,k_bot-1
167 tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
169 if (k_bot > k_top)
then
173 if ((kz < nz_in) .and. (kz > 1))
call &
174 find_limited_slope(tr_1d, z_edges, sl_tr, kz)
176 tr(i,j,k) = tr(i,j,k) + wt(kz) * &
177 (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
186 if ((e(k) > z_edges(1)) .and. (z_edges(nz_in+1) > e(k+1)))
then
187 tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
188 (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
189 (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
191 elseif (e(k) > z_edges(1))
then
192 tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
193 (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
195 elseif (z_edges(nz_in) > e(k+1))
then
196 tr(i,j,k) = ((e(k) - z_edges(nz_in+1)) * tr(i,j,k) + &
197 (z_edges(nz_in+1) - e(k+1)) * tr_1d(nz_in)) / &
203 do k=1,nz ; tr(i,j,k) = landval ;
enddo
209 do i=is,ie ; htot(i) = 0.0 ;
enddo
210 do k=1,nz ;
do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ;
enddo ;
enddo
212 do i=is,ie ;
if (g%mask2dT(i,j)*htot(i) > 0.0)
then
214 dilate = (g%bathyT(i,j) - 0.0) / htot(i)
215 e(nz+1) = -g%bathyT(i,j)
216 do k=nz,1,-1 ; e(k) = e(k+1) + dilate * h(i,j,k) ;
enddo
219 do k=1,nz_in ; tr_1d(k) = tr_in(i,j,k) ;
enddo
222 if (e(k+1) > z_edges(1))
then
224 elseif (z_edges(nz_in) > e(k))
then
225 tr(i,j,k) = tr_1d(nz_in)
227 call find_overlap(z_edges, e(k), e(k+1), nz_in-1, &
228 k_bot, k_top, k_bot, wt, z1, z2)
231 if (k_top < nz_in)
then
232 tr(i,j,k) = wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
233 (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
235 tr(i,j,k) = wt(kz)*tr_1d(nz_in)
237 do kz=k_top+1,k_bot-1
238 tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*(tr_1d(kz) + tr_1d(kz+1))
240 if (k_bot > k_top)
then
242 tr(i,j,k) = tr(i,j,k) + wt(kz)*0.5*((tr_1d(kz) + tr_1d(kz+1)) + &
243 (tr_1d(kz+1) - tr_1d(kz))*(z2(kz)+z1(kz)))
248 if ((e(k) > z_edges(1)) .and. (z_edges(nz_in) > e(k+1)))
then
249 tr(i,j,k) = (((e(k) - z_edges(1)) * tr_1d(1) + &
250 (z_edges(1) - z_edges(nz_in)) * tr(i,j,k)) + &
251 (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
253 elseif (e(k) > z_edges(1))
then
254 tr(i,j,k) = ((e(k) - z_edges(1)) * tr_1d(1) + &
255 (z_edges(1) - e(k+1)) * tr(i,j,k)) / &
257 elseif (z_edges(nz_in) > e(k+1))
then
258 tr(i,j,k) = ((e(k) - z_edges(nz_in)) * tr(i,j,k) + &
259 (z_edges(nz_in) - e(k+1)) * tr_1d(nz_in)) / &
265 do k=1,nz ; tr(i,j,k) = landval ;
enddo
270 deallocate(tr_in) ;
deallocate(tr_1d) ;
deallocate(z_edges)
271 deallocate(wt) ;
deallocate(z1) ;
deallocate(z2)
273 tracer_z_init = .true.