15 implicit none ;
private
17 #include <MOM_memory.h>
30 function tracer_z_init(tr, h, filename, tr_name, G, US, missing_val, land_val)
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)
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)
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 &
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 &
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)
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)
279 subroutine read_z_edges(filename, tr_name, z_edges, nz_out, has_edges, &
280 use_missing, missing, scale)
281 character(len=*),
intent(in) :: filename
282 character(len=*),
intent(in) :: tr_name
283 real,
dimension(:),
allocatable, &
284 intent(out) :: z_edges
285 integer,
intent(out) :: nz_out
286 logical,
intent(out) :: has_edges
288 logical,
intent(inout) :: use_missing
290 real,
intent(inout) :: missing
291 real,
intent(in) :: scale
295 character(len=32) :: mdl
296 character(len=120) :: dim_name, edge_name, tr_msg, dim_msg
298 integer :: ncid, status, intid, tr_id, layid, k
299 integer :: nz_edge, ndim, tr_dim_ids(NF90_MAX_VAR_DIMS)
301 mdl =
"MOM_tracer_Z_init read_Z_edges: "
302 tr_msg = trim(tr_name)//
" in "//trim(filename)
304 status = nf90_open(filename, nf90_nowrite, ncid)
305 if (status /= nf90_noerr)
then
306 call mom_error(warning,mdl//
" Difficulties opening "//trim(filename)//&
307 " - "//trim(nf90_strerror(status)))
311 status = nf90_inq_varid(ncid, tr_name, tr_id)
312 if (status /= nf90_noerr)
then
313 call mom_error(warning,mdl//
" Difficulties finding variable "//&
314 trim(tr_msg)//
" - "//trim(nf90_strerror(status)))
315 nz_out = -1 ; status = nf90_close(ncid) ;
return
317 status = nf90_inquire_variable(ncid, tr_id, ndims=ndim, dimids=tr_dim_ids)
318 if (status /= nf90_noerr)
then
319 call mom_error(warning,mdl//
" cannot inquire about "//trim(tr_msg))
320 elseif ((ndim < 3) .or. (ndim > 4))
then
321 call mom_error(warning,mdl//
" "//trim(tr_msg)//&
322 " has too many or too few dimensions.")
323 nz_out = -1 ; status = nf90_close(ncid) ;
return
326 if (.not.use_missing)
then
328 status = nf90_get_att(ncid, tr_id,
"missing_value", missing)
329 if (status /= nf90_noerr) use_missing = .true.
333 status = nf90_inquire_dimension(ncid, tr_dim_ids(3), dim_name, len=nz_out)
334 if (status /= nf90_noerr)
then
335 call mom_error(warning,mdl//
" cannot inquire about dimension(3) of "//&
339 dim_msg = trim(dim_name)//
" in "//trim(filename)
340 status = nf90_inq_varid(ncid, dim_name, layid)
341 if (status /= nf90_noerr)
then
342 call mom_error(warning,mdl//
" Difficulties finding variable "//&
343 trim(dim_msg)//
" - "//trim(nf90_strerror(status)))
344 nz_out = -1 ; status = nf90_close(ncid) ;
return
347 status = nf90_get_att(ncid, layid,
"edges", edge_name)
348 if (status /= nf90_noerr)
then
349 call mom_mesg(mdl//
" "//trim(dim_msg)//&
350 " has no readable edges attribute - "//trim(nf90_strerror(status)))
354 status = nf90_inq_varid(ncid, edge_name, intid)
355 if (status /= nf90_noerr)
then
356 call mom_error(warning,mdl//
" Difficulties finding edge variable "//&
357 trim(edge_name)//
" in "//trim(filename)//
" - "//trim(nf90_strerror(status)))
362 nz_edge = nz_out ;
if (has_edges) nz_edge = nz_out+1
363 allocate(z_edges(nz_edge)) ; z_edges(:) = 0.0
365 if (nz_out < 1)
return
369 dim_msg = trim(edge_name)//
" in "//trim(filename)
370 status = nf90_get_var(ncid, intid, z_edges)
371 if (status /= nf90_noerr)
then
372 call mom_error(warning,mdl//
" Difficulties reading variable "//&
373 trim(dim_msg)//
" - "//trim(nf90_strerror(status)))
374 nz_out = -1 ; status = nf90_close(ncid) ;
return
377 status = nf90_get_var(ncid, layid, z_edges)
378 if (status /= nf90_noerr)
then
379 call mom_error(warning,mdl//
" Difficulties reading variable "//&
380 trim(dim_msg)//
" - "//trim(nf90_strerror(status)))
381 nz_out = -1 ; status = nf90_close(ncid) ;
return
385 status = nf90_close(ncid)
386 if (status /= nf90_noerr)
call mom_error(warning, mdl// &
387 " Difficulties closing "//trim(filename)//
" - "//trim(nf90_strerror(status)))
391 if (z_edges(1) < z_edges(2))
then
392 do k=1,nz_edge ; z_edges(k) = -z_edges(k) ;
enddo
396 do k=2,nz_edge ;
if (z_edges(k) >= z_edges(k-1)) monotonic = .false. ;
enddo
397 if (.not.monotonic) &
398 call mom_error(warning,mdl//
" "//trim(dim_msg)//
" is not monotonic.")
400 if (scale /= 1.0)
then ;
do k=1,nz_edge ; z_edges(k) = scale*z_edges(k) ;
enddo ;
endif
415 subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
416 real,
dimension(:),
intent(in) :: e
417 real,
intent(in) :: Z_top
418 real,
intent(in) :: Z_bot
419 integer,
intent(in) :: k_max
420 integer,
intent(in) :: k_start
421 integer,
intent(inout) :: k_top
423 integer,
intent(inout) :: k_bot
425 real,
dimension(:),
intent(out) :: wt
426 real,
dimension(:),
intent(out) :: z1
429 real,
dimension(:),
intent(out) :: z2
433 real :: Ih, e_c, tot_wt, I_totwt
436 do k=k_start,k_max ;
if (e(k+1)<z_top)
exit ;
enddo
442 if (e(k+1)<=z_bot)
then
443 wt(k) = 1.0 ; k_bot = k
444 ih = 0.0 ;
if (e(k) /= e(k+1)) ih = 1.0 / (e(k)-e(k+1))
445 e_c = 0.5*(e(k)+e(k+1))
446 z1(k) = (e_c - min(e(k),z_top)) * ih
447 z2(k) = (e_c - z_bot) * ih
449 wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k)
450 if (e(k) /= e(k+1))
then
451 z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k), z_top)) / (e(k)-e(k+1))
452 else ; z1(k) = -0.5 ;
endif
456 if (e(k+1)<=z_bot)
then
458 wt(k) = e(k) - z_bot ; z1(k) = -0.5
459 if (e(k) /= e(k+1))
then
460 z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
461 else ; z2(k) = 0.5 ;
endif
463 wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
465 tot_wt = tot_wt + wt(k)
469 i_totwt = 1.0 / tot_wt
470 do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ;
enddo
479 real,
dimension(:),
intent(in) :: val
480 real,
dimension(:),
intent(in) :: e
481 real,
intent(out) :: slope
482 integer,
intent(in) :: k
486 d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
487 if (((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) .or. (d1*d2 <= 0.0))
then
490 slope = (d1**2*(val(k+1) - val(k)) + d2**2*(val(k) - val(k-1))) * &
491 ((e(k) - e(k+1)) / (d1*d2*(d1+d2)))
494 slope = sign(1.0,slope) * min(abs(slope), &
495 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
496 2.0*(val(k) - min(val(k-1),val(k),val(k+1))))