New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

trcsms_fabm.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC/FABM – NEMO

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90 @ 11979

Last change on this file since 11979 was 10728, checked in by dford, 5 years ago
Changes for assimilation of biogeochemical variables. See https://code.metoffice.gov.uk/trac/utils/ticket/174.
File size: 22.4 KB

Line
1	MODULE trcsms_fabm
2	!!======================================================================
3	!! * MODULE trcsms_fabm *
4	!! TOP : Main module of the FABM tracers
5	!!======================================================================
6	!! History : 1.0 ! 2015-04 (PML) Original code
7	!!----------------------------------------------------------------------
8	#if defined key_fabm
9	!!----------------------------------------------------------------------
10	!! 'key_fabm' FABM tracers
11	!!----------------------------------------------------------------------
12	!! trc_sms_fabm : FABM model main routine
13	!! trc_sms_fabm_alloc : allocate arrays specific to FABM sms
14	!!----------------------------------------------------------------------
15	USE par_trc ! TOP parameters
16	USE oce_trc ! Ocean variables
17	USE trc ! TOP variables
18	USE trcbc
19	USE trd_oce
20	USE trdtrc
21	#if defined key_trdtrc && defined key_iomput
22	USE trdtrc_oce
23	#endif
24
25	USE oce, only: tsn ! Needed?
26	USE sbc_oce, only: lk_oasis
27	USE dom_oce
28	USE zdf_oce
29	!USE iom
30	USE xios
31	USE cpl_oasis3
32	USE st2D_fabm
33	USE inputs_fabm
34	USE vertical_movement_fabm
35	USE zdfmxl
36	USE asmbgc, ONLY: mld_choice_bgc
37	USE lbclnk
38
39	!USE fldread ! time interpolation
40
41	USE fabm
42
43	IMPLICIT NONE
44
45	# include "domzgr_substitute.h90"
46	# include "vectopt_loop_substitute.h90"
47
48	PRIVATE
49
50	PUBLIC trc_sms_fabm ! called by trcsms.F90 module
51	PUBLIC trc_sms_fabm_alloc ! called by trcini_fabm.F90 module
52	PUBLIC trc_sms_fabm_check_mass
53	PUBLIC st2d_fabm_nxt ! 2D state intergration
54	PUBLIC compute_fabm ! Compute FABM sources, sinks and diagnostics
55
56	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: flux ! Cross-interface flux of pelagic variables (# m-2 s-1)
57
58	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ext ! Light extinction coefficient (m-1)
59
60	! Work array for mass aggregation
61	REAL(wp), ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: current_total
62
63
64	! Arrays for environmental variables
65	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: prn,rho
66	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: taubot
67
68	! repair counters
69	INTEGER :: repair_interior_count,repair_surface_count,repair_bottom_count
70
71	! state check type
72	TYPE type_state
73	LOGICAL :: valid
74	LOGICAL :: repaired
75	END TYPE
76
77	REAL(wp), PUBLIC :: daynumber_in_year
78
79	TYPE (type_bulk_variable_id),SAVE :: swr_id
80
81	!!----------------------------------------------------------------------
82	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
83	!! $Id$
84	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
85	!!----------------------------------------------------------------------
86	CONTAINS
87
88	SUBROUTINE trc_sms_fabm( kt )
89	!!----------------------------------------------------------------------
90	!! * trc_sms_fabm *
91	!!
92	!! ** Purpose : main routine of FABM model
93	!!
94	!! ** Method : -
95	!!----------------------------------------------------------------------
96	!
97	INTEGER, INTENT(in) :: kt ! ocean time-step index
98	INTEGER :: jn
99	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrfabm
100
101	!!----------------------------------------------------------------------
102	!
103	IF( nn_timing == 1 ) CALL timing_start('trc_sms_fabm')
104	!
105	IF(lwp) WRITE(numout,*)
106	IF(lwp) WRITE(numout,'(a,i0,a,i4.4,a,i2.2,a,i2.2,a,i5,a)') &
107	' trc_sms_fabm: FABM model, iteration ',kt,' ', &
108	nyear,'-',nmonth,'-',nday,' ',nsec_day," secs"
109	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
110
111	CALL update_inputs( kt )
112
113	CALL compute_fabm
114
115	CALL compute_vertical_movement( kt )
116
117	CALL st2d_fabm_nxt( kt )
118
119	CALL asmdiags_fabm( kt )
120
121	IF( l_trdtrc ) CALL wrk_alloc( jpi, jpj, jpk, ztrfabm )
122
123	CALL trc_bc_read ( kt ) ! tracers: surface and lateral Boundary Conditions
124	CALL trc_rnf_fabm ( kt ) ! River forcings
125
126	IF( l_trdtrc ) THEN ! Save the trends in the mixed layer
127	DO jn = jp_fabm0, jp_fabm1
128	ztrfabm(:,:,:) = tra(:,:,:,jn)
129	CALL trd_trc( ztrfabm, jn, jptra_sms, kt ) ! save trends
130	END DO
131	CALL wrk_dealloc( jpi, jpj, jpk, ztrfabm )
132	END IF
133	!
134	IF( nn_timing == 1 ) CALL timing_stop('trc_sms_fabm')
135
136	END SUBROUTINE trc_sms_fabm
137
138	SUBROUTINE asmdiags_fabm( kt )
139	INTEGER, INTENT(IN) :: kt
140	INTEGER :: ji,jj,jk,jkmax
141	REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld
142
143	IF (kt == nittrc000) THEN
144	MLD_MAX(:,:) = 0.0
145	ENDIF
146	PGROW_AVG(:,:) = 0.0
147	PLOSS_AVG(:,:) = 0.0
148	PHYT_AVG(:,:) = 0.0
149
150	pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow)
151	ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss)
152
153	SELECT CASE( mld_choice_bgc )
154	CASE ( 1 ) ! Turbocline/mixing depth [W points]
155	zmld(:,:) = hmld(:,:)
156	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
157	zmld(:,:) = hmlp(:,:)
158	CASE ( 3 ) ! Kara MLD [Interpolated]
159	#if defined key_karaml
160	IF ( ln_kara ) THEN
161	zmld(:,:) = hmld_kara(:,:)
162	ELSE
163	CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
164	& ' but ln_kara=.false.' )
165	ENDIF
166	#else
167	CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
168	& ' but is not defined' )
169	#endif
170	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
171	zmld(:,:) = hmld_tref(:,:)
172	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
173	zmld(:,:) = hmlpt(:,:)
174	END SELECT
175
176	DO jj = 2, jpjm1
177	DO ji = 2, jpim1
178	!
179	jkmax = jpk-1
180	DO jk = jpk-1, 1, -1
181	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
182	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
183	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
184	jkmax = jk
185	ENDIF
186	END DO
187	!
188	DO jk = 1, jkmax
189	PHYT_AVG(ji,jj) = PHYT_AVG(ji,jj) + &
190	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
191	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
192	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
193	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
194	IF ( pgrow_3d(ji,jj,jk) .GT. 0.0 ) THEN
195	PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) + &
196	& pgrow_3d(ji,jj,jk)
197	ENDIF
198	IF ( ploss_3d(ji,jj,jk) .GT. 0.0 ) THEN
199	PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) + &
200	& ploss_3d(ji,jj,jk)
201	ENDIF
202	END DO
203
204	PHYT_AVG(ji,jj) = PHYT_AVG(ji,jj) / REAL(jkmax)
205	PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) / REAL(jkmax)
206	PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) / REAL(jkmax)
207
208	IF ( zmld(ji,jj) .GT. MLD_MAX(ji,jj) ) THEN
209	MLD_MAX(ji,jj) = zmld(ji,jj)
210	ENDIF
211	!
212	END DO
213	END DO
214
215	PHYT_AVG(:,:) = PHYT_AVG(:,:) * tmask(:,:,1)
216	PGROW_AVG(:,:) = PGROW_AVG(:,:) * tmask(:,:,1)
217	PLOSS_AVG(:,:) = PLOSS_AVG(:,:) * tmask(:,:,1)
218	MLD_MAX(:,:) = MLD_MAX(:,:) * tmask(:,:,1)
219
220	END SUBROUTINE asmdiags_fabm
221
222	SUBROUTINE compute_fabm()
223	INTEGER :: ji,jj,jk,jn
224	TYPE(type_state) :: valid_state
225	REAL(wp) :: zalfg,zztmpx,zztmpy
226
227	! Validate current model state (setting argument to .TRUE. enables repair=clipping)
228	valid_state = check_state(.TRUE.)
229	IF (.NOT.valid_state%valid) THEN
230	WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!"
231	#if defined key_iomput
232	CALL xios_finalize ! end mpp communications with xios
233	IF( lk_oasis ) CALL cpl_finalize ! end coupling and mpp communications with OASIS
234	#else
235	IF( lk_oasis ) THEN
236	CALL cpl_finalize ! end coupling and mpp communications with OASIS
237	ELSE
238	IF( lk_mpp ) CALL mppstop ! end mpp communications
239	ENDIF
240	#endif
241	END IF
242	IF (valid_state%repaired) THEN
243	WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count
244	WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count
245	WRITE(numout,*) "Total bottom repairs up to now on process",narea,":",repair_bottom_count
246	ENDIF
247
248	! Compute the now hydrostatic pressure
249	! copied from istate.F90
250	! ------------------------------------
251
252	zalfg = 0.5e-4 * grav ! FABM wants dbar, convert from Pa
253
254	rho = rau0 * ( 1. + rhd )
255
256	prn(:,:,1) = 10.1325 + zalfg * fse3t(:,:,1) * rho(:,:,1)
257
258	daynumber_in_year=(fjulday-fjulstartyear+1)*1._wp
259
260	DO jk = 2, jpk ! Vertical integration from the surface
261	prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( &
262	fse3t(:,:,jk-1) * rho(:,:,jk-1) &
263	+ fse3t(:,:,jk) * rho(:,:,jk) )
264	END DO
265
266	! Bottom stress
267	taubot(:,:) = 0._wp
268	DO jj = 2, jpjm1
269	DO ji = fs_2, fs_jpim1 ! vector opt.
270	zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) &
271	& + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) )
272	zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) &
273	& + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) )
274	taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
275	!
276	ENDDO
277	ENDDO
278	! Compute light extinction
279	DO jk=1,jpk
280	DO jj=1,jpj
281	call fabm_get_light_extinction(model,1,jpi,jj,jk,ext)
282	END DO
283	END DO
284
285	! Compute light field (stored among FABM's internal diagnostics)
286	DO jj=1,jpj
287	DO ji=1,jpi
288	call fabm_get_light(model,1,jpk,ji,jj)
289	END DO
290	END DO
291
292	! TODO: retrieve 3D shortwave and store in etot3
293
294	! Zero rate array of interface-attached state variables
295	fabm_st2Da = 0._wp
296
297	! Compute interfacial source terms and fluxes
298	DO jj=1,jpj
299	! Process bottom (fabm_do_bottom increments rather than sets, so zero flux array first)
300	flux = 0._wp
301	CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:))
302	DO jn=1,jp_fabm
303	DO ji=1,jpi
304	! Divide bottom fluxes by height of bottom layer and add to source terms.
305	! TODO: is there perhaps an existing variable for fse3t(ji,jj,mbkt(ji,jj))??
306	tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) = tra(ji,jj,mbkt(ji,jj),jp_fabm_m1+jn) + flux(ji,jn)/fse3t(ji,jj,mbkt(ji,jj))
307	END DO
308	END DO
309
310	! Process surface (fabm_do_surface increments rather than sets, so zero flux array first)
311	flux = 0._wp
312	CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface))
313	DO jn=1,jp_fabm
314	! Divide surface fluxes by height of surface layer and add to source terms.
315	tra(:,jj,1,jp_fabm_m1+jn) = tra(:,jj,1,jp_fabm_m1+jn) + flux(:,jn)/fse3t(:,jj,1)
316	END DO
317	END DO
318
319	! Compute interior source terms (NB fabm_do increments rather than sets)
320	DO jk=1,jpk
321	DO jj=1,jpj
322	CALL fabm_do(model,1,jpi,jj,jk,tra(:,jj,jk,jp_fabm0:jp_fabm1))
323	END DO
324	END DO
325	END SUBROUTINE compute_fabm
326
327	FUNCTION check_state(repair) RESULT(exit_state)
328	LOGICAL, INTENT(IN) :: repair
329	TYPE(type_state) :: exit_state
330
331	INTEGER :: jj,jk
332	LOGICAL :: valid_int,valid_sf,valid_bt
333
334	exit_state%valid = .TRUE.
335	exit_state%repaired =.FALSE.
336	DO jk=1,jpk
337	DO jj=1,jpj
338	CALL fabm_check_state(model,1,jpi,jj,jk,repair,valid_int)
339	IF (repair.AND..NOT.valid_int) THEN
340	repair_interior_count = repair_interior_count + 1
341	exit_state%repaired = .TRUE.
342	END IF
343	IF (.NOT.(valid_int.OR.repair)) exit_state%valid = .FALSE.
344	END DO
345	END DO
346	DO jj=1,jpj
347	CALL fabm_check_surface_state(model,1,jpi,jj,repair,valid_sf)
348	IF (repair.AND..NOT.valid_sf) THEN
349	repair_surface_count = repair_surface_count + 1
350	exit_state%repaired = .TRUE.
351	END IF
352	IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
353	CALL fabm_check_bottom_state(model,1,jpi,jj,repair,valid_bt)
354	IF (repair.AND..NOT.valid_bt) THEN
355	repair_bottom_count = repair_bottom_count + 1
356	exit_state%repaired = .TRUE.
357	END IF
358	IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
359	END DO
360	END FUNCTION
361
362	SUBROUTINE trc_sms_fabm_check_mass()
363	REAL(wp) :: total(SIZE(model%conserved_quantities))
364	INTEGER :: jk,jj,jn
365
366	total = 0._wp
367
368	DO jk=1,jpk
369	DO jj=1,jpj
370	CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total)
371	DO jn=1,SIZE(model%conserved_quantities)
372	total(jn) = total(jn) + SUM(cvol(:,jj,jk)current_total(:,jn)tmask_i(:,jj))
373	END DO
374	END DO
375	END DO
376
377	DO jj=1,jpj
378	CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total)
379	DO jn=1,SIZE(model%conserved_quantities)
380	total(jn) = total(jn) + SUM(e1e2t(:,jj)current_total(:,jn)tmask_i(:,jj))
381	END DO
382	END DO
383
384	IF( lk_mpp ) CALL mpp_sum(total,SIZE(model%conserved_quantities))
385
386	DO jn=1,SIZE(model%conserved_quantities)
387	IF(lwp) WRITE(numout,) 'FABM '//TRIM(model%conserved_quantities(jn)%name),total(jn),TRIM(model%conserved_quantities(jn)%units)//'m3'
388	END DO
389
390	END SUBROUTINE trc_sms_fabm_check_mass
391
392	SUBROUTINE st2d_fabm_nxt( kt )
393	!!----------------------------------------------------------------------
394	!! * st2d_fabm_nxt *
395	!!
396	!! ** Purpose : routine to integrate 2d states in time
397	!!
398	!! ** Method : based on integration of 3D passive tracer fields
399	!! implemented in TOP_SRC/TRP/trcnxt.F90, plus
400	!! tra_nxt_fix in OPA_SRC/TRA/tranxt.F90. Similar to
401	!! time integration of sea surface height in
402	!! OPA_SRC/DYN/sshwzv.F90.
403	!!----------------------------------------------------------------------
404	!
405	INTEGER, INTENT(in) :: kt ! ocean time-step index
406	REAL(wp) :: z2dt
407	INTEGER :: jn
408
409	!!----------------------------------------------------------------------
410	!
411	IF( neuler == 0 .AND. kt == nittrc000 ) THEN
412	z2dt = rdt ! set time step size (Euler)
413	ELSE
414	z2dt = 2._wp * rdt ! set time step size (Leapfrog)
415	ENDIF
416
417	! Forward Euler time step to compute "now"
418	DO jn=1,jp_fabm_surface+jp_fabm_bottom
419	fabm_st2Da(:,:,jn) = (fabm_st2db(:,:,jn) + z2dt * fabm_st2da(:,:,jn)) * tmask(:,:,1)
420	ENDDO
421
422	IF( neuler == 0 .AND. kt == nittrc000 ) THEN ! Euler time-stepping at first time-step
423	! ! (only swap)
424	fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
425	!
426	ELSE
427	! Update now state + Asselin filter time stepping
428	fabm_st2Db(:,:,:) = (1._wp - 2._wpatfp) fabm_st2Dn(:,:,:) + &
429	atfp * ( fabm_st2Db(:,:,:) + fabm_st2Da(:,:,:) )
430	fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
431	ENDIF
432
433	END SUBROUTINE st2d_fabm_nxt
434
435	INTEGER FUNCTION trc_sms_fabm_alloc()
436	INTEGER :: jj,jk,jn
437	!!----------------------------------------------------------------------
438	!! * ROUTINE trc_sms_fabm_alloc *
439	!!----------------------------------------------------------------------
440	!
441	! ALLOCATE here the arrays specific to FABM
442	ALLOCATE( lk_rad_fabm(jp_fabm))
443	ALLOCATE( prn(jpi, jpj, jpk))
444	ALLOCATE( rho(jpi, jpj, jpk))
445	ALLOCATE( taubot(jpi, jpj))
446	! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc )
447
448	! Allocate arrays to hold state for surface-attached and bottom-attached state variables
449	ALLOCATE(fabm_st2Dn(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
450	ALLOCATE(fabm_st2Da(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
451	ALLOCATE(fabm_st2Db(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
452
453	! Work array to hold surface and bottom fluxes
454	ALLOCATE(flux(jpi,jp_fabm))
455
456	! Work array to hold extinction coefficients
457	ALLOCATE(ext(jpi))
458	ext=0._wp
459
460	! Allocate work arrays for vertical movement
461	ALLOCATE(w_ct(jpi,jpk,jp_fabm))
462	ALLOCATE(w_if(jpk,jp_fabm))
463	ALLOCATE(zwgt_if(jpk,jp_fabm))
464	ALLOCATE(flux_if(jpk,jp_fabm))
465	ALLOCATE(current_total(jpi,SIZE(model%conserved_quantities)))
466	#if defined key_trdtrc && defined key_iomput
467	IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm))
468	IF( lk_trdtrc ) ALLOCATE(tr_inp(jpi,jpj,jpk))
469	#endif
470
471	trc_sms_fabm_alloc = 0 ! set to zero if no array to be allocated
472	!
473	IF( trc_sms_fabm_alloc /= 0 ) CALL ctl_warn('trc_sms_fabm_alloc : failed to allocate arrays')
474	!
475
476	! Make FABM aware of diagnostics that are not needed [not included in output]
477	DO jn=1,size(model%diagnostic_variables)
478	!model%diagnostic_variables(jn)%save = iom_use(model%diagnostic_variables(jn)%name)
479	END DO
480	DO jn=1,size(model%horizontal_diagnostic_variables)
481	!model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name)
482	END DO
483
484	! Provide FABM with domain extents [after this, the save attribute of diagnostic variables can no longe change!]
485	call fabm_set_domain(model,jpi, jpj, jpk)
486
487	! Provide FABM with the vertical indices of the surface and bottom, and the land-sea mask.
488	call model%set_bottom_index(mbkt) ! NB mbkt extents should match dimension lengths provided to fabm_set_domain
489	call model%set_surface_index(1)
490	call fabm_set_mask(model,tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to fabm_set_domain
491
492	! Send pointers to state data to FABM
493	do jn=1,jp_fabm
494	call fabm_link_bulk_state_data(model,jn,trn(:,:,:,jp_fabm_m1+jn))
495	end do
496	DO jn=1,jp_fabm_surface
497	CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn))
498	END DO
499	DO jn=1,jp_fabm_bottom
500	CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn))
501	END DO
502
503	! Send pointers to environmental data to FABM
504	call fabm_link_bulk_data(model,standard_variables%temperature,tsn(:,:,:,jp_tem))
505	call fabm_link_bulk_data(model,standard_variables%practical_salinity,tsn(:,:,:,jp_sal))
506	call fabm_link_bulk_data(model,standard_variables%density,rho(:,:,:))
507	call fabm_link_bulk_data(model,standard_variables%pressure,prn)
508	call fabm_link_horizontal_data(model,standard_variables%bottom_stress,taubot(:,:))
509	! correct target for cell thickness depends on NEMO configuration:
510	#ifdef key_vvl
511	call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_n)
512	#else
513	call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_0)
514	#endif
515	call fabm_link_horizontal_data(model,standard_variables%latitude,gphit)
516	call fabm_link_horizontal_data(model,standard_variables%longitude,glamt)
517	call fabm_link_scalar_data(model,standard_variables%number_of_days_since_start_of_the_year,daynumber_in_year)
518	call fabm_link_horizontal_data(model,standard_variables%wind_speed,wndm(:,:))
519	call fabm_link_horizontal_data(model,standard_variables%surface_downwelling_shortwave_flux,qsr(:,:))
520	call fabm_link_horizontal_data(model,standard_variables%bottom_depth_below_geoid,bathy(:,:))
521
522	swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux)
523
524	! Obtain user-specified input variables (read from NetCDF file)
525	call link_inputs
526	call update_inputs( nit000, .false. )
527
528	! Check whether FABM has all required data
529	call fabm_check_ready(model)
530
531	! Initialize state
532	DO jj=1,jpj
533	CALL fabm_initialize_surface_state(model,1,jpi,jj)
534	CALL fabm_initialize_bottom_state(model,1,jpi,jj)
535	END DO
536	DO jk=1,jpk
537	DO jj=1,jpj
538	CALL fabm_initialize_state(model,1,jpi,jj,jk)
539	END DO
540	END DO
541
542	! Set mask for negativity corrections to the relevant states
543	lk_rad_fabm = .FALSE.
544	DO jn=1,jp_fabm
545	IF (model%state_variables(jn)%minimum.ge.0) THEN
546	lk_rad_fabm(jn)=.TRUE.
547	IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%state_variables(jn)%name)//' activated.'
548	END IF
549	END DO
550
551	! Mask land points in states with zeros, not nice, but coherent
552	! with NEMO "convention":
553	DO jn=jp_fabm0,jp_fabm1
554	WHERE (tmask==0._wp)
555	trn(:,:,:,jn)=0._wp
556	END WHERE
557	END DO
558	DO jn=1,jp_fabm_surface+jp_fabm_bottom
559	WHERE (tmask(:,:,1)==0._wp)
560	fabm_st2Dn(:,:,jn)=0._wp
561	END WHERE
562	END DO
563
564	! Copy initial condition for interface-attached state variables to "previous" state field
565	! NB NEMO does this itself for pelagic state variables (trb) in TOP_SRC/trcini.F90.
566	fabm_st2Db = fabm_st2Dn
567
568	! Initialise repair counters
569	repair_interior_count = 0
570	repair_surface_count = 0
571	repair_bottom_count = 0
572
573	END FUNCTION trc_sms_fabm_alloc
574
575	#else
576	!!----------------------------------------------------------------------
577	!! Dummy module No FABM model
578	!!----------------------------------------------------------------------
579	CONTAINS
580	SUBROUTINE trc_sms_fabm( kt ) ! Empty routine
581	INTEGER, INTENT( in ) :: kt
582	WRITE(,) 'trc_sms_fabm: You should not have seen this print! error?', kt
583	END SUBROUTINE trc_sms_fabm
584	#endif
585
586	!!======================================================================
587	END MODULE trcsms_fabm

Note: See TracBrowser for help on using the repository browser.

Download in other formats: