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/CO6_KD490_amm7_oper_fabm_hem08/NEMOGCM/NEMO/TOP_SRC/FABM – NEMO

source: branches/UKMO/CO6_KD490_amm7_oper_fabm_hem08/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90 @ 9336

Last change on this file since 9336 was 9336, checked in by dford, 6 years ago

Add growth and loss diagnostics.

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