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, 20 months ago

Changes for assimilation of biogeochemical variables. See https://code.metoffice.gov.uk/trac/utils/ticket/174.

File size: 22.4 KB
Line 
1MODULE 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   !!----------------------------------------------------------------------
86CONTAINS
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._wp*atfp) * 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   !!----------------------------------------------------------------------
579CONTAINS
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   !!======================================================================
587END MODULE trcsms_fabm
Note: See TracBrowser for help on using the repository browser.