Line data Source code
1 : #ifdef ncdf
2 : #define USE_NETCDF
3 : #endif
4 : !=======================================================================
5 : !
6 : ! Reads and interpolates forcing data for biogeochemistry
7 : !
8 : ! authors: Nicole Jeffery, LANL
9 : ! Elizabeth C. Hunke, LANL
10 : !
11 : module ice_forcing_bgc
12 :
13 : use ice_kinds_mod
14 : use ice_blocks, only: nx_block, ny_block
15 : use ice_domain_size, only: max_blocks
16 : use ice_communicate, only: my_task, master_task
17 : use ice_calendar, only: dt, istep, msec, mday, mmonth
18 : use ice_fileunits, only: nu_diag
19 : use ice_arrays_column, only: restore_bgc, &
20 : bgc_data_dir, fe_data_type
21 : use ice_constants, only: c0, p1
22 : use ice_constants, only: field_loc_center, field_type_scalar
23 : use ice_exit, only: abort_ice
24 : use ice_forcing, only: bgc_data_type
25 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
26 : use icepack_intfc, only: icepack_nspint_3bd, icepack_max_aero, &
27 : icepack_max_algae, icepack_max_doc, icepack_max_dic
28 : use icepack_intfc, only: icepack_query_tracer_flags, &
29 : icepack_query_parameters, icepack_query_parameters, & ! LCOV_EXCL_LINE
30 : icepack_query_tracer_indices
31 :
32 : implicit none
33 : private
34 : public :: get_forcing_bgc, get_atm_bgc, fzaero_data, alloc_forcing_bgc, &
35 : init_bgc_data, faero_data, faero_default, fiso_default
36 :
37 : integer (kind=int_kind) :: &
38 : bgcrecnum = 0 ! old record number (save between steps)
39 :
40 : real (kind=dbl_kind), dimension(:,:,:), allocatable, public :: &
41 : nitdat , & ! data value toward which nitrate is restored ! LCOV_EXCL_LINE
42 : sildat ! data value toward which silicate is restored
43 :
44 : real (kind=dbl_kind), dimension(:,:,:,:), allocatable, public :: &
45 : nit_data, & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
46 : sil_data
47 :
48 : !=======================================================================
49 :
50 : contains
51 :
52 : !=======================================================================
53 : !
54 : ! Allocate space for forcing_bgc variables
55 : !
56 0 : subroutine alloc_forcing_bgc
57 :
58 : integer (int_kind) :: ierr
59 :
60 : allocate( &
61 : nitdat (nx_block,ny_block,max_blocks), & ! data value toward which nitrate is restored ! LCOV_EXCL_LINE
62 : sildat (nx_block,ny_block,max_blocks), & ! data value toward which silicate is restored ! LCOV_EXCL_LINE
63 : nit_data(nx_block,ny_block,2,max_blocks), & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
64 : sil_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
65 0 : stat=ierr)
66 0 : if (ierr/=0) call abort_ice('(alloc_forcing_bgc): Out of memory')
67 :
68 0 : end subroutine alloc_forcing_bgc
69 :
70 : !=======================================================================
71 : !
72 : ! Read and interpolate annual climatologies of silicate and nitrate.
73 : ! Restore model quantities to data if desired.
74 : !
75 : ! author: Elizabeth C. Hunke, LANL
76 :
77 0 : subroutine get_forcing_bgc
78 :
79 : use ice_blocks, only: block, get_block
80 : use ice_domain, only: nblocks, blocks_ice
81 : use ice_arrays_column, only: ocean_bio_all
82 : use ice_calendar, only: yday
83 : ! use ice_flux, only: sss
84 : use ice_flux_bgc, only: sil, nit
85 : use ice_forcing, only: trestore, trest, fyear, &
86 : read_clim_data_nc, interpolate_data, & ! LCOV_EXCL_LINE
87 : interp_coeff_monthly, interp_coeff, & ! LCOV_EXCL_LINE
88 : read_data_nc_point, c1intp, c2intp
89 :
90 : integer (kind=int_kind) :: &
91 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
92 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
93 : ixm,ixp, ixx , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
94 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
95 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
96 : midmonth , & ! middle day of month ! LCOV_EXCL_LINE
97 : recnum , & ! record number ! LCOV_EXCL_LINE
98 : dataloc , & ! = 1 for data located in middle of time interval ! LCOV_EXCL_LINE
99 : ! = 2 for date located at end of time interval
100 : ks ! bgc tracer index (bio_index_o)
101 :
102 : character (char_len_long) :: &
103 : met_file, & ! netcdf filename ! LCOV_EXCL_LINE
104 : fieldname ! field name in netcdf file
105 :
106 : real (kind=dbl_kind), dimension(2), save :: &
107 : sil_data_p , & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
108 : nit_data_p ! field values at 2 temporal data points
109 :
110 : real (kind=dbl_kind) :: &
111 : secday , & ! number of seconds in day ! LCOV_EXCL_LINE
112 0 : sec1hr ! number of seconds in 1 hour
113 :
114 : logical (kind=log_kind) :: readm, read1, tr_bgc_Nit, tr_bgc_Sil
115 :
116 : character (char_len_long) :: & ! input data file names
117 : nit_file , & ! nitrate input file ! LCOV_EXCL_LINE
118 : sil_file ! silicate input file
119 :
120 : type (block) :: &
121 : this_block ! block information for current block
122 :
123 : character(len=*), parameter :: subname = '(get_forcing_bgc)'
124 :
125 0 : call icepack_query_parameters(secday_out=secday)
126 : call icepack_query_tracer_flags(tr_bgc_Nit_out=tr_bgc_Nit, &
127 0 : tr_bgc_Sil_out=tr_bgc_Sil)
128 0 : call icepack_warnings_flush(nu_diag)
129 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
130 0 : file=__FILE__, line=__LINE__)
131 :
132 0 : if (trim(bgc_data_type) == 'clim') then
133 :
134 0 : nit_file = 'nitrate_climatologyWOA_gx1v6f_20150107.nc'
135 : !'nitrate_WOA2005_surface_monthly' ! gx1 only
136 0 : sil_file = 'silicate_climatologyWOA_gx1v6f_20150107.nc'
137 : !'silicate_WOA2005_surface_monthly' ! gx1 only
138 :
139 0 : nit_file = trim(bgc_data_dir)//'/'//trim(nit_file)
140 0 : sil_file = trim(bgc_data_dir)//'/'//trim(sil_file)
141 :
142 0 : if (my_task == master_task .and. istep == 1) then
143 0 : if (tr_bgc_Sil) then
144 0 : write (nu_diag,*) ' '
145 0 : write (nu_diag,*) 'silicate data interpolated to timestep:'
146 0 : write (nu_diag,*) trim(sil_file)
147 : endif
148 0 : if (tr_bgc_Nit) then
149 0 : write (nu_diag,*) ' '
150 0 : write (nu_diag,*) 'nitrate data interpolated to timestep:'
151 0 : write (nu_diag,*) trim(nit_file)
152 0 : if (restore_bgc) write (nu_diag,*) &
153 0 : 'bgc restoring timescale (days) =', trestore
154 : endif
155 : endif ! my_task, istep
156 :
157 : !-------------------------------------------------------------------
158 : ! monthly data
159 : !
160 : ! Assume that monthly data values are located in the middle of the
161 : ! month.
162 : !-------------------------------------------------------------------
163 :
164 0 : midmonth = 15 ! data is given on 15th of every month
165 : !!! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
166 :
167 : ! Compute record numbers for surrounding months
168 0 : maxrec = 12
169 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
170 0 : ixp = mod(mmonth, maxrec) + 1
171 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
172 0 : if (mday < midmonth) ixp = -99
173 :
174 : ! Determine whether interpolation will use values 1:2 or 2:3
175 : ! recslot = 2 means we use values 1:2, with the current value (2)
176 : ! in the second slot
177 : ! recslot = 1 means we use values 2:3, with the current value (2)
178 : ! in the first slot
179 0 : recslot = 1 ! latter half of month
180 0 : if (mday < midmonth) recslot = 2 ! first half of month
181 :
182 : ! Find interpolation coefficients
183 0 : call interp_coeff_monthly (recslot)
184 :
185 0 : readm = .false.
186 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
187 :
188 : endif ! 'clim prep'
189 :
190 : !-------------------------------------------------------------------
191 : ! Read two monthly silicate values and interpolate.
192 : ! Restore toward interpolated value.
193 : !-------------------------------------------------------------------
194 :
195 0 : if (trim(bgc_data_type)=='clim' .AND. tr_bgc_Sil) then
196 : ! call read_clim_data (readm, 0, ixm, mmonth, ixp, &
197 : ! sil_file, sil_data, & ! LCOV_EXCL_LINE
198 : ! field_loc_center, field_type_scalar)
199 0 : fieldname = 'silicate'
200 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
201 : sil_file, fieldname, sil_data, & ! LCOV_EXCL_LINE
202 0 : field_loc_center, field_type_scalar)
203 0 : call interpolate_data (sil_data, sildat)
204 :
205 0 : if (istep == 1 .or. .NOT. restore_bgc) then
206 :
207 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
208 0 : do iblk = 1, nblocks
209 :
210 0 : this_block = get_block(blocks_ice(iblk),iblk)
211 0 : ilo = this_block%ilo
212 0 : ihi = this_block%ihi
213 0 : jlo = this_block%jlo
214 0 : jhi = this_block%jhi
215 :
216 0 : do j = jlo, jhi
217 0 : do i = ilo, ihi
218 :
219 0 : sil(i,j,iblk) = sildat(i,j,iblk)
220 0 : ks = 2*icepack_max_algae + icepack_max_doc + 3 + icepack_max_dic
221 0 : ocean_bio_all(i,j,ks,iblk) = sil(i,j,iblk) !Sil
222 : enddo
223 : enddo
224 : enddo
225 :
226 : !$OMP END PARALLEL DO
227 0 : elseif (restore_bgc) then
228 :
229 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
230 0 : do iblk = 1, nblocks
231 :
232 0 : this_block = get_block(blocks_ice(iblk),iblk)
233 0 : ilo = this_block%ilo
234 0 : ihi = this_block%ihi
235 0 : jlo = this_block%jlo
236 0 : jhi = this_block%jhi
237 :
238 0 : do j = jlo, jhi
239 0 : do i = ilo, ihi
240 :
241 : sil(i,j,iblk) = sil(i,j,iblk) &
242 0 : + (sildat(i,j,iblk)-sil(i,j,iblk))*dt/trest
243 0 : ks = 2*icepack_max_algae + icepack_max_doc + 3 + icepack_max_dic
244 0 : ocean_bio_all(i,j,ks,iblk) = sil(i,j,iblk) !Sil
245 : enddo
246 : enddo
247 : enddo
248 : !$OMP END PARALLEL DO
249 : endif !restore
250 0 : elseif (tr_bgc_Sil) then ! bgc_data_type /= 'clim'
251 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
252 0 : do iblk = 1, nblocks
253 :
254 0 : this_block = get_block(blocks_ice(iblk),iblk)
255 0 : ilo = this_block%ilo
256 0 : ihi = this_block%ihi
257 0 : jlo = this_block%jlo
258 0 : jhi = this_block%jhi
259 :
260 0 : do j = jlo, jhi
261 0 : do i = ilo, ihi
262 :
263 0 : sil(i,j,iblk) = 25.0_dbl_kind
264 0 : ks = 2*icepack_max_algae + icepack_max_doc + 3 + icepack_max_dic
265 0 : ocean_bio_all(i,j,ks,iblk) = sil(i,j,iblk) !Sil
266 : enddo
267 : enddo
268 : enddo
269 : !$OMP END PARALLEL DO
270 :
271 : endif !tr_bgc_Sil
272 : !-------------------------------------------------------------------
273 : ! Read two monthly nitrate values and interpolate.
274 : ! Restore toward interpolated value.
275 : !-------------------------------------------------------------------
276 :
277 0 : if (trim(bgc_data_type)=='clim' .AND. tr_bgc_Nit) then
278 : ! call read_clim_data (readm, 0, ixm, mmonth, ixp, &
279 : ! nit_file, nit_data, & ! LCOV_EXCL_LINE
280 : ! field_loc_center, field_type_scalar)
281 0 : fieldname = 'nitrate'
282 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
283 : nit_file, fieldname, nit_data, & ! LCOV_EXCL_LINE
284 0 : field_loc_center, field_type_scalar)
285 0 : call interpolate_data (nit_data, nitdat)
286 :
287 0 : if (istep == 1 .or. .NOT. restore_bgc) then
288 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
289 0 : do iblk = 1, nblocks
290 :
291 0 : this_block = get_block(blocks_ice(iblk),iblk)
292 0 : ilo = this_block%ilo
293 0 : ihi = this_block%ihi
294 0 : jlo = this_block%jlo
295 0 : jhi = this_block%jhi
296 :
297 0 : do j = jlo, jhi
298 0 : do i = ilo, ihi
299 :
300 0 : nit(i,j,iblk) = nitdat(i,j,iblk)
301 0 : ks = icepack_max_algae + 1
302 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !nit
303 0 : ks = 2*icepack_max_algae + icepack_max_doc + 7 + icepack_max_dic
304 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !PON
305 : enddo
306 : enddo
307 : enddo
308 : !$OMP END PARALLEL DO
309 0 : elseif (restore_bgc ) then
310 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
311 0 : do iblk = 1, nblocks
312 :
313 0 : this_block = get_block(blocks_ice(iblk),iblk)
314 0 : ilo = this_block%ilo
315 0 : ihi = this_block%ihi
316 0 : jlo = this_block%jlo
317 0 : jhi = this_block%jhi
318 :
319 0 : do j = jlo, jhi
320 0 : do i = ilo, ihi
321 :
322 : nit(i,j,iblk) = nit(i,j,iblk) &
323 0 : + (nitdat(i,j,iblk)-nit(i,j,iblk))*dt/trest
324 0 : ks = icepack_max_algae + 1
325 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !nit
326 0 : ks = 2*icepack_max_algae + icepack_max_doc + 7 + icepack_max_dic
327 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !PON
328 : enddo
329 : enddo
330 : enddo
331 : !$OMP END PARALLEL DO
332 : endif !restore_bgc
333 :
334 : ! elseif (trim(nit_data_type) == 'sss' .AND. tr_bgc_Nit) then
335 : ! !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
336 : ! do iblk = 1, nblocks
337 :
338 : ! this_block = get_block(blocks_ice(iblk),iblk)
339 : ! ilo = this_block%ilo
340 : ! ihi = this_block%ihi
341 : ! jlo = this_block%jlo
342 : ! jhi = this_block%jhi
343 :
344 : ! do j = jlo, jhi
345 : ! do i = ilo, ihi
346 :
347 : ! nit(i,j,iblk) = sss(i,j,iblk)
348 : ! ks = icepack_max_algae + 1
349 : ! ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !nit
350 : ! ks = 2*icepack_max_algae + icepack_max_doc + 7 + icepack_max_dic
351 : ! ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !PON
352 : ! enddo
353 : ! enddo
354 : ! enddo
355 : ! !$OMP END PARALLEL DO
356 :
357 0 : elseif (tr_bgc_Nit) then ! bgc_data_type /= 'clim'
358 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
359 0 : do iblk = 1, nblocks
360 :
361 0 : this_block = get_block(blocks_ice(iblk),iblk)
362 0 : ilo = this_block%ilo
363 0 : ihi = this_block%ihi
364 0 : jlo = this_block%jlo
365 0 : jhi = this_block%jhi
366 :
367 0 : do j = jlo, jhi
368 0 : do i = ilo, ihi
369 :
370 0 : nit(i,j,iblk) = 12.0_dbl_kind
371 0 : ks = icepack_max_algae + 1
372 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !nit
373 0 : ks = 2*icepack_max_algae + icepack_max_doc + 7 + icepack_max_dic
374 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !PON
375 : enddo
376 : enddo
377 : enddo
378 : !$OMP END PARALLEL DO
379 :
380 : endif !tr_bgc_Nit
381 :
382 : !-------------------------------------------------------------------
383 : ! Data from Papdimitrious et al., 2007, Limnol. Oceanogr.
384 : ! and WOA at 68oS, 304.5oE :
385 : ! daily data located at the end of the 24-hour period.
386 : !-------------------------------------------------------------------
387 :
388 0 : if (trim(bgc_data_type) == 'ISPOL') then
389 :
390 0 : nit_file = trim(bgc_data_dir)//'nutrients_daily_ISPOL_WOA_field3.nc'
391 0 : sil_file = trim(bgc_data_dir)//'nutrients_daily_ISPOL_WOA_field3.nc'
392 :
393 0 : if (my_task == master_task .and. istep == 1) then
394 0 : if (tr_bgc_Sil) then
395 0 : write (nu_diag,*) ' '
396 0 : write (nu_diag,*) 'silicate data interpolated to timestep:'
397 0 : write (nu_diag,*) trim(sil_file)
398 : endif
399 0 : if (tr_bgc_Nit) then
400 0 : write (nu_diag,*) ' '
401 0 : write (nu_diag,*) 'nitrate data interpolated to timestep:'
402 0 : write (nu_diag,*) trim(nit_file)
403 0 : if (restore_bgc) write (nu_diag,*) &
404 0 : 'bgc restoring timescale (days) =', trestore
405 : endif
406 : endif ! my_task, istep
407 :
408 0 : dataloc = 2 ! data located at end of interval
409 0 : sec1hr = secday ! seconds in day
410 0 : maxrec = 365 !
411 :
412 : ! current record number
413 0 : recnum = int(yday)
414 :
415 : ! Compute record numbers for surrounding data (2 on each side)
416 0 : ixm = mod(recnum+maxrec-2,maxrec) + 1
417 0 : ixx = mod(recnum-1, maxrec) + 1
418 :
419 0 : recslot = 2
420 0 : ixp = -99
421 0 : call interp_coeff (recnum, recslot, sec1hr, dataloc)
422 :
423 0 : read1 = .false.
424 0 : if (istep==1 .or. bgcrecnum .ne. recnum) read1 = .true.
425 :
426 :
427 0 : if (tr_bgc_Sil) then
428 0 : met_file = sil_file
429 0 : fieldname= 'silicate'
430 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
431 : maxrec, met_file, fieldname, sil_data_p, & ! LCOV_EXCL_LINE
432 0 : field_loc_center, field_type_scalar)
433 :
434 0 : sil(:,:,:) = c1intp * sil_data_p(1) &
435 0 : + c2intp * sil_data_p(2)
436 : endif
437 :
438 0 : if (tr_bgc_Nit) then
439 0 : met_file = nit_file
440 0 : fieldname= 'nitrate'
441 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
442 : maxrec, met_file, fieldname, nit_data_p, & ! LCOV_EXCL_LINE
443 0 : field_loc_center, field_type_scalar)
444 :
445 0 : nit(:,:,:) = c1intp * nit_data_p(1) &
446 0 : + c2intp * nit_data_p(2)
447 : endif
448 :
449 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
450 0 : do iblk = 1, nblocks
451 :
452 0 : this_block = get_block(blocks_ice(iblk),iblk)
453 0 : ilo = this_block%ilo
454 0 : ihi = this_block%ihi
455 0 : jlo = this_block%jlo
456 0 : jhi = this_block%jhi
457 :
458 0 : do j = jlo, jhi
459 0 : do i = ilo, ihi
460 :
461 0 : ks = 2*icepack_max_algae + icepack_max_doc + 3 + icepack_max_dic
462 0 : ocean_bio_all(i,j,ks,iblk) = sil(i,j,iblk) !Sil
463 0 : ks = icepack_max_algae + 1
464 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !nit
465 0 : ks = 2*icepack_max_algae + icepack_max_doc + 7 + icepack_max_dic
466 0 : ocean_bio_all(i,j,ks,iblk) = nit(i,j,iblk) !PON
467 : enddo
468 : enddo
469 : enddo
470 : !$OMP END PARALLEL DO
471 :
472 : ! Save record number for next time step
473 0 : bgcrecnum = recnum
474 : endif
475 :
476 0 : end subroutine get_forcing_bgc
477 :
478 : !=======================================================================
479 : !
480 : ! author: Nicole Jeffery, LANL
481 :
482 0 : subroutine get_atm_bgc
483 :
484 : use ice_blocks, only: block, get_block
485 : use ice_domain, only: nblocks, blocks_ice
486 : use ice_domain_size, only: n_zaero
487 : use ice_flux_bgc, only: flux_bio_atm, faero_atm
488 :
489 : ! local variables
490 :
491 : integer (kind=int_kind) :: &
492 : i, j, nn , & ! horizontal indices ! LCOV_EXCL_LINE
493 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
494 : iblk ! block index
495 :
496 : logical (kind=log_kind) :: &
497 : tr_zaero
498 :
499 : integer (kind=int_kind), dimension(icepack_max_aero) :: &
500 : nlt_zaero
501 :
502 : type (block) :: &
503 : this_block ! block information for current block
504 :
505 : character(len=*), parameter :: subname = '(get_atm_bgc)'
506 :
507 : !-----------------------------------------------------------------
508 : ! initialize
509 : !-----------------------------------------------------------------
510 :
511 0 : call icepack_query_tracer_flags(tr_zaero_out=tr_zaero)
512 0 : call icepack_query_tracer_indices(nlt_zaero_out=nlt_zaero)
513 0 : call icepack_warnings_flush(nu_diag)
514 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
515 0 : file=__FILE__, line=__LINE__)
516 :
517 0 : flux_bio_atm(:,:,:,:) = c0
518 0 : if (tr_zaero) then
519 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,nn)
520 0 : do iblk = 1, nblocks
521 :
522 0 : this_block = get_block(blocks_ice(iblk),iblk)
523 0 : ilo = this_block%ilo
524 0 : ihi = this_block%ihi
525 0 : jlo = this_block%jlo
526 0 : jhi = this_block%jhi
527 :
528 0 : do nn = 1, n_zaero
529 0 : do j = jlo, jhi
530 0 : do i = ilo, ihi
531 0 : flux_bio_atm(i,j,nlt_zaero(nn),iblk) = faero_atm(i,j,nn,iblk)
532 : enddo
533 : enddo
534 : enddo
535 :
536 : enddo ! iblk
537 : !$OMP END PARALLEL DO
538 : endif
539 :
540 0 : end subroutine get_atm_bgc
541 :
542 : !=======================================================================
543 :
544 : ! constant values for atmospheric water isotopes
545 : !
546 : ! authors: David Bailey, NCAR
547 :
548 0 : subroutine fiso_default
549 :
550 : use ice_flux_bgc, only: fiso_atm
551 : character(len=*), parameter :: subname='(fiso_default)'
552 :
553 0 : fiso_atm(:,:,:,:) = 1.e-14_dbl_kind ! kg/m^2 s
554 :
555 0 : end subroutine fiso_default
556 :
557 : !=======================================================================
558 :
559 : ! constant values for atmospheric aerosols
560 : !
561 : ! authors: Elizabeth Hunke, LANL
562 :
563 0 : subroutine faero_default
564 :
565 : use ice_flux_bgc, only: faero_atm
566 :
567 : character(len=*), parameter :: subname = '(faero_default)'
568 :
569 0 : faero_atm(:,:,1,:) = 1.e-12_dbl_kind ! kg/m^2 s
570 0 : faero_atm(:,:,2,:) = 1.e-13_dbl_kind
571 0 : faero_atm(:,:,3,:) = 1.e-14_dbl_kind
572 0 : faero_atm(:,:,4,:) = 1.e-14_dbl_kind
573 0 : faero_atm(:,:,5,:) = 1.e-14_dbl_kind
574 0 : faero_atm(:,:,6,:) = 1.e-14_dbl_kind
575 :
576 0 : end subroutine faero_default
577 :
578 : !=======================================================================
579 :
580 : ! read atmospheric aerosols
581 : !
582 : ! authors: Elizabeth Hunke, LANL
583 :
584 0 : subroutine faero_data
585 :
586 : use ice_calendar, only: mmonth, mday, istep, msec
587 : use ice_domain_size, only: max_blocks
588 : use ice_blocks, only: nx_block, ny_block
589 : use ice_flux_bgc, only: faero_atm
590 : use ice_forcing, only: interp_coeff_monthly, read_clim_data_nc, interpolate_data
591 :
592 : ! local parameters
593 :
594 : real (kind=dbl_kind), dimension(:,:,:,:), allocatable, &
595 : save :: & ! LCOV_EXCL_LINE
596 : aero1_data , & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
597 : aero2_data , & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
598 : aero3_data ! field values at 2 temporal data points
599 :
600 : character (char_len_long) :: &
601 : aero_file, & ! netcdf filename ! LCOV_EXCL_LINE
602 : fieldname ! field name in netcdf file
603 :
604 : integer (kind=int_kind) :: &
605 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
606 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
607 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
608 : midmonth ! middle day of month
609 :
610 : logical (kind=log_kind) :: readm
611 :
612 : character(len=*), parameter :: subname = '(faero_data)'
613 :
614 0 : allocate( aero1_data(nx_block,ny_block,2,max_blocks), &
615 : aero2_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
616 0 : aero3_data(nx_block,ny_block,2,max_blocks) )
617 :
618 :
619 : !-------------------------------------------------------------------
620 : ! monthly data
621 : !
622 : ! Assume that monthly data values are located in the middle of the
623 : ! month.
624 : !-------------------------------------------------------------------
625 :
626 0 : midmonth = 15 ! data is given on 15th of every month
627 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
628 :
629 : ! Compute record numbers for surrounding months
630 0 : maxrec = 12
631 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
632 0 : ixp = mod(mmonth, maxrec) + 1
633 0 : if (mday >= midmonth) ixm = 99 ! other two points will be used
634 0 : if (mday < midmonth) ixp = 99
635 :
636 : ! Determine whether interpolation will use values 1:2 or 2:3
637 : ! recslot = 2 means we use values 1:2, with the current value (2)
638 : ! in the second slot
639 : ! recslot = 1 means we use values 2:3, with the current value (2)
640 : ! in the first slot
641 0 : recslot = 1 ! latter half of month
642 0 : if (mday < midmonth) recslot = 2 ! first half of month
643 :
644 : ! Find interpolation coefficients
645 0 : call interp_coeff_monthly (recslot)
646 :
647 : ! Read 2 monthly values
648 0 : readm = .false.
649 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
650 :
651 : ! aero_file = trim(atm_data_dir)//'faero.nc'
652 0 : aero_file = '/usr/projects/climate/eclare/DATA/gx1v3/faero.nc'
653 :
654 0 : fieldname='faero_atm001'
655 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
656 : aero_file, fieldname, aero1_data, & ! LCOV_EXCL_LINE
657 0 : field_loc_center, field_type_scalar)
658 :
659 0 : fieldname='faero_atm002'
660 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
661 : aero_file, fieldname, aero2_data, & ! LCOV_EXCL_LINE
662 0 : field_loc_center, field_type_scalar)
663 :
664 0 : fieldname='faero_atm003'
665 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
666 : aero_file, fieldname, aero3_data, & ! LCOV_EXCL_LINE
667 0 : field_loc_center, field_type_scalar)
668 :
669 0 : call interpolate_data (aero1_data, faero_atm(:,:,1,:)) ! W/m^2 s
670 0 : call interpolate_data (aero2_data, faero_atm(:,:,2,:))
671 0 : call interpolate_data (aero3_data, faero_atm(:,:,3,:))
672 :
673 0 : where (faero_atm(:,:,:,:) > 1.e20) faero_atm(:,:,:,:) = c0
674 :
675 0 : deallocate( aero1_data, aero2_data, aero3_data )
676 :
677 0 : end subroutine faero_data
678 :
679 : !=======================================================================
680 :
681 : ! read atmospheric aerosols
682 : !
683 : ! authors: Elizabeth Hunke, LANL
684 :
685 0 : subroutine fzaero_data
686 :
687 : use ice_blocks, only: nx_block, ny_block
688 : use ice_flux_bgc, only: faero_atm
689 : use ice_forcing, only: interp_coeff_monthly, read_clim_data_nc, interpolate_data
690 :
691 : ! local parameters
692 :
693 : real (kind=dbl_kind), dimension(:,:,:,:), allocatable, &
694 : save :: & ! LCOV_EXCL_LINE
695 : aero_data ! field values at 2 temporal data points
696 :
697 : character (char_len_long) :: &
698 : aero_file, & ! netcdf filename ! LCOV_EXCL_LINE
699 : fieldname ! field name in netcdf file
700 :
701 : integer (kind=int_kind) :: &
702 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
703 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
704 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
705 : midmonth ! middle day of month
706 :
707 : integer (kind=int_kind), dimension(icepack_max_aero) :: &
708 : nlt_zaero
709 :
710 : logical (kind=log_kind) :: readm
711 :
712 : character(len=*), parameter :: subname = '(fzaero_data)'
713 :
714 0 : call icepack_query_tracer_indices(nlt_zaero_out=nlt_zaero)
715 0 : call icepack_warnings_flush(nu_diag)
716 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
717 0 : file=__FILE__, line=__LINE__)
718 :
719 0 : allocate( aero_data(nx_block,ny_block,2,max_blocks) )
720 :
721 : !-------------------------------------------------------------------
722 : ! monthly data
723 : !
724 : ! Assume that monthly data values are located in the middle of the
725 : ! month.
726 : !-------------------------------------------------------------------
727 :
728 0 : midmonth = 15 ! data is given on 15th of every month
729 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
730 :
731 : ! Compute record numbers for surrounding months
732 0 : maxrec = 12
733 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
734 0 : ixp = mod(mmonth, maxrec) + 1
735 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
736 0 : if (mday < midmonth) ixp = -99
737 :
738 : ! Determine whether interpolation will use values 1:2 or 2:3
739 : ! recslot = 2 means we use values 1:2, with the current value (2)
740 : ! in the second slot
741 : ! recslot = 1 means we use values 2:3, with the current value (2)
742 : ! in the first slot
743 0 : recslot = 1 ! latter half of month
744 0 : if (mday < midmonth) recslot = 2 ! first half of month
745 :
746 : ! Find interpolation coefficients
747 0 : call interp_coeff_monthly (recslot)
748 :
749 : ! Read 2 monthly values
750 0 : readm = .false.
751 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
752 :
753 : ! aero_file = trim(atm_data_dir)//'faero.nc'
754 : ! Cam5 monthly total black carbon deposition on the gx1 grid"
755 0 : aero_file = '/usr/projects/climate/njeffery/DATA/CAM/Hailong_Wang/Cam5_bc_monthly_popgrid.nc'
756 :
757 0 : fieldname='bcd'
758 : call read_clim_data_nc (readm, 0, ixm, mmonth, ixp, &
759 : aero_file, fieldname, aero_data, & ! LCOV_EXCL_LINE
760 0 : field_loc_center, field_type_scalar)
761 :
762 :
763 0 : call interpolate_data (aero_data, faero_atm(:,:,nlt_zaero(1),:)) ! kg/m^2/s
764 :
765 0 : where (faero_atm(:,:,nlt_zaero(1),:) > 1.e20) faero_atm(:,:,nlt_zaero(1),:) = c0
766 :
767 0 : deallocate( aero_data )
768 :
769 0 : end subroutine fzaero_data
770 :
771 : !=======================================================================
772 :
773 : ! Initialize ocean iron from file
774 : !
775 : ! authors: Nicole Jeffery, LANL
776 :
777 0 : subroutine init_bgc_data (fed1,fep1)
778 :
779 : use ice_read_write, only: ice_open_nc, ice_read_nc, ice_close_nc
780 :
781 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(inout) :: &
782 : fed1, & ! first dissolved iron pool (nM) ! LCOV_EXCL_LINE
783 : fep1 ! first particulate iron pool (nM)
784 :
785 : ! local parameters
786 :
787 : integer (kind=int_kind) :: &
788 : fid ! file id for netCDF file
789 :
790 : logical (kind=log_kind) :: diag
791 :
792 : character (char_len_long) :: &
793 : iron_file, & ! netcdf filename ! LCOV_EXCL_LINE
794 : fieldname ! field name in netcdf file
795 :
796 : character(len=*), parameter :: subname = '(init_bgc_data)'
797 :
798 : !-------------------------------------------------------------------
799 : ! Annual average data from Tagliabue, 2012 (top 50 m average
800 : ! poisson grid filled on gx1v6
801 : !-------------------------------------------------------------------
802 :
803 0 : if (trim(fe_data_type) == 'clim') then
804 0 : diag = .true. ! write diagnostic information
805 0 : iron_file = trim(bgc_data_dir)//'dFe_50m_annual_Tagliabue_gx1.nc'
806 :
807 0 : if (my_task == master_task) then
808 0 : write (nu_diag,*) ' '
809 0 : write (nu_diag,*) 'Dissolved iron ocean concentrations from:'
810 0 : write (nu_diag,*) trim(iron_file)
811 0 : call ice_open_nc(iron_file,fid)
812 : endif
813 :
814 0 : fieldname='dFe'
815 : ! Currently only first fed value is read
816 0 : call ice_read_nc(fid,1,fieldname,fed1,diag)
817 0 : where ( fed1(:,:,:) > 1.e20) fed1(:,:,:) = p1
818 :
819 0 : if (my_task == master_task) call ice_close_nc(fid)
820 :
821 0 : diag = .true. ! write diagnostic information
822 0 : iron_file = trim(bgc_data_dir)//'pFe_bathy_gx1.nc'
823 :
824 0 : if (my_task == master_task) then
825 0 : write (nu_diag,*) ' '
826 0 : write (nu_diag,*) 'Particulate iron ocean concentrations from:'
827 0 : write (nu_diag,*) trim(iron_file)
828 0 : call ice_open_nc(iron_file,fid)
829 : endif
830 :
831 0 : fieldname='pFe'
832 : ! Currently only first fep value is read
833 0 : call ice_read_nc(fid,1,fieldname,fep1,diag)
834 0 : where ( fep1(:,:,:) > 1.e20) fep1(:,:,:) = p1
835 :
836 0 : if (my_task == master_task) call ice_close_nc(fid)
837 :
838 : endif
839 :
840 0 : end subroutine init_bgc_data
841 :
842 : !=======================================================================
843 :
844 : end module ice_forcing_bgc
845 :
846 : !=======================================================================
|
