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

source: branches/UKMO/AMM15_v3_6_STABLE_package_bgc_updates/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90 @ 10241

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

Add visibility diagnostic (can't currently be done via fabm.yaml without changing the FABM code).
File size: 19.2 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
36   !USE fldread         !  time interpolation
37
38   USE fabm
39   USE par_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      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrfabm )
118
119      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions
120      CALL trc_rnf_fabm ( kt ) ! River forcings
121     
122      visib(:,:,:) = 1.7 / fabm_get_bulk_diagnostic_data(model, jp_fabm_xeps)
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 compute_fabm()
137      INTEGER :: ji,jj,jk,jn
138      TYPE(type_state) :: valid_state
139      REAL(wp) :: zalfg,zztmpx,zztmpy
140
141      ! Validate current model state (setting argument to .TRUE. enables repair=clipping)
142      valid_state = check_state(.TRUE.)
143      IF (.NOT.valid_state%valid) THEN
144         WRITE(numout,*) "Invalid value in FABM encountered in area ",narea,"!!!"
145#if defined key_iomput
146         CALL xios_finalize                ! end mpp communications with xios
147         IF( lk_oasis ) CALL cpl_finalize    ! end coupling and mpp communications with OASIS
148#else
149         IF( lk_oasis ) THEN
150            CALL cpl_finalize              ! end coupling and mpp communications with OASIS
151         ELSE
152            IF( lk_mpp )   CALL mppstop    ! end mpp communications
153         ENDIF
154#endif
155      END IF
156      IF (valid_state%repaired) THEN
157         WRITE(numout,*) "Total interior repairs up to now on process",narea,":",repair_interior_count
158         WRITE(numout,*) "Total surface repairs up to now on process",narea,":",repair_surface_count
159         WRITE(numout,*) "Total bottom repairs up to now on process",narea,":",repair_bottom_count
160      ENDIF
161
162      ! Compute the now hydrostatic pressure
163      ! copied from istate.F90
164      ! ------------------------------------
165
166      zalfg = 0.5e-4 * grav ! FABM wants dbar, convert from Pa
167
168      rho = rau0 * ( 1. + rhd )
169
170      prn(:,:,1) = 10.1325 + zalfg * fse3t(:,:,1) * rho(:,:,1)
171
172      daynumber_in_year=(fjulday-fjulstartyear+1)*1._wp
173
174      DO jk = 2, jpk                                              ! Vertical integration from the surface
175         prn(:,:,jk) = prn(:,:,jk-1) + zalfg * ( &
176                     fse3t(:,:,jk-1) * rho(:,:,jk-1)  &
177                     + fse3t(:,:,jk) * rho(:,:,jk) )
178      END DO
179
180      ! Bottom stress
181      taubot(:,:) = 0._wp
182      DO jj = 2, jpjm1
183         DO ji = fs_2, fs_jpim1   ! vector opt.
184               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
185                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )
186               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
187                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )
188               taubot(ji,jj) = 0.5_wp * rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
189               !
190         ENDDO
191      ENDDO
192      ! Compute light extinction
193      DO jk=1,jpk
194          DO jj=1,jpj
195            call fabm_get_light_extinction(model,1,jpi,jj,jk,ext)
196         END DO
197      END DO
198
199      ! Compute light field (stored among FABM's internal diagnostics)
200      DO jj=1,jpj
201          DO ji=1,jpi
202            call fabm_get_light(model,1,jpk,ji,jj)
203         END DO
204      END DO
205
206      ! TODO: retrieve 3D shortwave and store in etot3
207
208      ! Zero rate array of interface-attached state variables
209      fabm_st2Da = 0._wp
210
211      ! Compute interfacial source terms and fluxes
212      DO jj=1,jpj
213         ! Process bottom (fabm_do_bottom increments rather than sets, so zero flux array first)
214         flux = 0._wp
215         CALL fabm_do_bottom(model,1,jpi,jj,flux,fabm_st2Da(:,jj,jp_fabm_surface+1:))
216         DO jn=1,jp_fabm
217             DO ji=1,jpi
218                 ! Divide bottom fluxes by height of bottom layer and add to source terms.
219                 ! TODO: is there perhaps an existing variable for fse3t(ji,jj,mbkt(ji,jj))??
220                 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))
221             END DO
222         END DO
223
224         ! Process surface (fabm_do_surface increments rather than sets, so zero flux array first)
225         flux = 0._wp
226         CALL fabm_do_surface(model,1,jpi,jj,flux,fabm_st2Da(:,jj,1:jp_fabm_surface))
227         DO jn=1,jp_fabm
228             ! Divide surface fluxes by height of surface layer and add to source terms.
229             tra(:,jj,1,jp_fabm_m1+jn) = tra(:,jj,1,jp_fabm_m1+jn) + flux(:,jn)/fse3t(:,jj,1)
230         END DO
231      END DO
232
233      ! Compute interior source terms (NB fabm_do increments rather than sets)
234      DO jk=1,jpk
235          DO jj=1,jpj
236              CALL fabm_do(model,1,jpi,jj,jk,tra(:,jj,jk,jp_fabm0:jp_fabm1))
237          END DO
238      END DO
239   END SUBROUTINE compute_fabm
240
241   FUNCTION check_state(repair) RESULT(exit_state)
242      LOGICAL, INTENT(IN) :: repair
243      TYPE(type_state) :: exit_state
244
245      INTEGER             :: jj,jk
246      LOGICAL             :: valid_int,valid_sf,valid_bt
247
248      exit_state%valid = .TRUE.
249      exit_state%repaired =.FALSE.
250      DO jk=1,jpk
251         DO jj=1,jpj
252            CALL fabm_check_state(model,1,jpi,jj,jk,repair,valid_int)
253            IF (repair.AND..NOT.valid_int) THEN
254               repair_interior_count = repair_interior_count + 1
255               exit_state%repaired = .TRUE.
256            END IF
257            IF (.NOT.(valid_int.OR.repair)) exit_state%valid = .FALSE.
258         END DO
259      END DO
260      DO jj=1,jpj
261         CALL fabm_check_surface_state(model,1,jpi,jj,repair,valid_sf)
262         IF (repair.AND..NOT.valid_sf) THEN
263            repair_surface_count = repair_surface_count + 1
264            exit_state%repaired = .TRUE.
265         END IF
266         IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
267         CALL fabm_check_bottom_state(model,1,jpi,jj,repair,valid_bt)
268         IF (repair.AND..NOT.valid_bt) THEN
269            repair_bottom_count = repair_bottom_count + 1
270            exit_state%repaired = .TRUE.
271         END IF
272         IF (.NOT.(valid_sf.AND.valid_bt).AND..NOT.repair) exit_state%valid = .FALSE.
273      END DO
274   END FUNCTION
275
276   SUBROUTINE trc_sms_fabm_check_mass()
277      REAL(wp) :: total(SIZE(model%conserved_quantities))
278      INTEGER :: jk,jj,jn
279
280      total = 0._wp
281
282      DO jk=1,jpk
283         DO jj=1,jpj
284            CALL fabm_get_conserved_quantities(model,1,jpi,jj,jk,current_total)
285            DO jn=1,SIZE(model%conserved_quantities)
286               total(jn) = total(jn) + SUM(cvol(:,jj,jk)*current_total(:,jn)*tmask_i(:,jj))
287            END DO
288         END DO
289      END DO
290
291      DO jj=1,jpj
292         CALL fabm_get_horizontal_conserved_quantities(model,1,jpi,jj,current_total)
293         DO jn=1,SIZE(model%conserved_quantities)
294            total(jn) = total(jn) + SUM(e1e2t(:,jj)*current_total(:,jn)*tmask_i(:,jj))
295         END DO
296      END DO
297
298      IF( lk_mpp ) CALL mpp_sum(total,SIZE(model%conserved_quantities))
299
300      DO jn=1,SIZE(model%conserved_quantities)
301         IF(lwp) WRITE(numout,*) 'FABM '//TRIM(model%conserved_quantities(jn)%name),total(jn),TRIM(model%conserved_quantities(jn)%units)//'*m3'
302      END DO
303
304   END SUBROUTINE trc_sms_fabm_check_mass
305
306   SUBROUTINE st2d_fabm_nxt( kt )
307      !!----------------------------------------------------------------------
308      !!                     ***  st2d_fabm_nxt  ***
309      !!
310      !! ** Purpose :   routine to integrate 2d states in time
311      !!
312      !! ** Method  :   based on integration of 3D passive tracer fields
313      !!                implemented in TOP_SRC/TRP/trcnxt.F90, plus
314      !!                tra_nxt_fix in OPA_SRC/TRA/tranxt.F90. Similar to
315      !!                time integration of sea surface height in
316      !!                OPA_SRC/DYN/sshwzv.F90.
317      !!----------------------------------------------------------------------
318      !
319      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
320      REAL(wp) :: z2dt
321      INTEGER :: jn
322
323!!----------------------------------------------------------------------
324      !
325      IF( neuler == 0 .AND. kt == nittrc000 ) THEN
326          z2dt = rdt                  ! set time step size (Euler)
327      ELSE
328          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
329      ENDIF
330
331      ! Forward Euler time step to compute "now"
332      DO jn=1,jp_fabm_surface+jp_fabm_bottom
333         fabm_st2Da(:,:,jn) = (fabm_st2db(:,:,jn) + z2dt * fabm_st2da(:,:,jn)) * tmask(:,:,1)
334      ENDDO
335
336      IF( neuler == 0 .AND. kt == nittrc000 ) THEN        ! Euler time-stepping at first time-step
337         !                                                ! (only swap)
338         fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
339         !
340      ELSE
341         ! Update now state + Asselin filter time stepping
342         fabm_st2Db(:,:,:) = (1._wp - 2._wp*atfp) * fabm_st2Dn(:,:,:) + &
343             atfp * ( fabm_st2Db(:,:,:) + fabm_st2Da(:,:,:) )
344         fabm_st2Dn(:,:,:) = fabm_st2Da(:,:,:)
345      ENDIF
346
347   END SUBROUTINE st2d_fabm_nxt
348
349   INTEGER FUNCTION trc_sms_fabm_alloc()
350      INTEGER :: jj,jk,jn
351      !!----------------------------------------------------------------------
352      !!              ***  ROUTINE trc_sms_fabm_alloc  ***
353      !!----------------------------------------------------------------------
354      !
355      ! ALLOCATE here the arrays specific to FABM
356      ALLOCATE( lk_rad_fabm(jp_fabm))
357      ALLOCATE( prn(jpi, jpj, jpk))
358      ALLOCATE( rho(jpi, jpj, jpk))
359      ALLOCATE( taubot(jpi, jpj))
360      ! ALLOCATE( tab(...) , STAT=trc_sms_fabm_alloc )
361
362      ! Allocate arrays to hold state for surface-attached and bottom-attached state variables
363      ALLOCATE(fabm_st2Dn(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
364      ALLOCATE(fabm_st2Da(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
365      ALLOCATE(fabm_st2Db(jpi, jpj, jp_fabm_surface+jp_fabm_bottom))
366
367      ! Work array to hold surface and bottom fluxes
368      ALLOCATE(flux(jpi,jp_fabm))
369
370      ! Work array to hold extinction coefficients
371      ALLOCATE(ext(jpi))
372      ext=0._wp
373
374      ! Allocate work arrays for vertical movement
375      ALLOCATE(w_ct(jpi,jpk,jp_fabm))
376      ALLOCATE(w_if(jpk,jp_fabm))
377      ALLOCATE(zwgt_if(jpk,jp_fabm))
378      ALLOCATE(flux_if(jpk,jp_fabm))
379      ALLOCATE(current_total(jpi,SIZE(model%conserved_quantities)))
380#if defined key_trdtrc && defined key_iomput
381      IF( lk_trdtrc ) ALLOCATE(tr_vmv(jpi,jpj,jpk,jp_fabm))
382      IF( lk_trdtrc ) ALLOCATE(tr_inp(jpi,jpj,jpk))
383#endif
384
385      trc_sms_fabm_alloc = 0      ! set to zero if no array to be allocated
386      !
387      IF( trc_sms_fabm_alloc /= 0 ) CALL ctl_warn('trc_sms_fabm_alloc : failed to allocate arrays')
388      !
389
390      ! Make FABM aware of diagnostics that are not needed [not included in output]
391      DO jn=1,size(model%diagnostic_variables)
392          !model%diagnostic_variables(jn)%save = iom_use(model%diagnostic_variables(jn)%name)
393      END DO
394      DO jn=1,size(model%horizontal_diagnostic_variables)
395          !model%horizontal_diagnostic_variables(jn)%save = iom_use(model%horizontal_diagnostic_variables(jn)%name)
396      END DO
397
398      ! Provide FABM with domain extents [after this, the save attribute of diagnostic variables can no longe change!]
399      call fabm_set_domain(model,jpi, jpj, jpk)
400
401      ! Provide FABM with the vertical indices of the surface and bottom, and the land-sea mask.
402      call model%set_bottom_index(mbkt)  ! NB mbkt extents should match dimension lengths provided to fabm_set_domain
403      call model%set_surface_index(1)
404      call fabm_set_mask(model,tmask,tmask(:,:,1)) ! NB tmask extents should match dimension lengths provided to fabm_set_domain
405
406      ! Send pointers to state data to FABM
407      do jn=1,jp_fabm
408         call fabm_link_bulk_state_data(model,jn,trn(:,:,:,jp_fabm_m1+jn))
409      end do
410      DO jn=1,jp_fabm_surface
411         CALL fabm_link_surface_state_data(model,jn,fabm_st2Dn(:,:,jn))
412      END DO
413      DO jn=1,jp_fabm_bottom
414         CALL fabm_link_bottom_state_data(model,jn,fabm_st2Dn(:,:,jp_fabm_surface+jn))
415      END DO
416
417      ! Send pointers to environmental data to FABM
418      call fabm_link_bulk_data(model,standard_variables%temperature,tsn(:,:,:,jp_tem))
419      call fabm_link_bulk_data(model,standard_variables%practical_salinity,tsn(:,:,:,jp_sal))
420      call fabm_link_bulk_data(model,standard_variables%density,rho(:,:,:))
421      call fabm_link_bulk_data(model,standard_variables%pressure,prn)
422      call fabm_link_horizontal_data(model,standard_variables%bottom_stress,taubot(:,:))
423      ! correct target for cell thickness depends on NEMO configuration:
424#ifdef key_vvl
425      call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_n)
426#else
427      call fabm_link_bulk_data(model,standard_variables%cell_thickness,e3t_0)
428#endif
429      call fabm_link_horizontal_data(model,standard_variables%latitude,gphit)
430      call fabm_link_horizontal_data(model,standard_variables%longitude,glamt)
431      call fabm_link_scalar_data(model,standard_variables%number_of_days_since_start_of_the_year,daynumber_in_year)
432      call fabm_link_horizontal_data(model,standard_variables%wind_speed,wndm(:,:))
433      call fabm_link_horizontal_data(model,standard_variables%surface_downwelling_shortwave_flux,qsr(:,:))
434      call fabm_link_horizontal_data(model,standard_variables%bottom_depth_below_geoid,bathy(:,:))
435
436      swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux)
437
438      ! Obtain user-specified input variables (read from NetCDF file)
439      call link_inputs
440      call update_inputs( nit000, .false. )
441
442      ! Check whether FABM has all required data
443      call fabm_check_ready(model)
444
445      ! Initialize state
446      DO jj=1,jpj
447         CALL fabm_initialize_surface_state(model,1,jpi,jj)
448         CALL fabm_initialize_bottom_state(model,1,jpi,jj)
449      END DO
450      DO jk=1,jpk
451         DO jj=1,jpj
452            CALL fabm_initialize_state(model,1,jpi,jj,jk)
453         END DO
454      END DO
455
456      ! Set mask for negativity corrections to the relevant states
457      lk_rad_fabm = .FALSE.
458      DO jn=1,jp_fabm
459        IF (model%state_variables(jn)%minimum.ge.0) THEN
460          lk_rad_fabm(jn)=.TRUE.
461          IF(lwp) WRITE(numout,*) 'FABM clipping for '//TRIM(model%state_variables(jn)%name)//' activated.'
462        END IF
463      END DO
464
465      ! Mask land points in states with zeros, not nice, but coherent
466      ! with NEMO "convention":
467      DO jn=jp_fabm0,jp_fabm1
468        WHERE (tmask==0._wp)
469          trn(:,:,:,jn)=0._wp
470        END WHERE
471      END DO
472      DO jn=1,jp_fabm_surface+jp_fabm_bottom
473        WHERE (tmask(:,:,1)==0._wp)
474          fabm_st2Dn(:,:,jn)=0._wp
475        END WHERE
476      END DO
477
478      ! Copy initial condition for interface-attached state variables to "previous" state field
479      ! NB NEMO does this itself for pelagic state variables (trb) in TOP_SRC/trcini.F90.
480      fabm_st2Db = fabm_st2Dn
481
482      ! Initialise repair counters
483      repair_interior_count = 0
484      repair_surface_count = 0
485      repair_bottom_count = 0
486
487   END FUNCTION trc_sms_fabm_alloc
488
489#else
490   !!----------------------------------------------------------------------
491   !!   Dummy module                                        No FABM model
492   !!----------------------------------------------------------------------
493CONTAINS
494   SUBROUTINE trc_sms_fabm( kt )             ! Empty routine
495      INTEGER, INTENT( in ) ::   kt
496      WRITE(*,*) 'trc_sms_fabm: You should not have seen this print! error?', kt
497   END SUBROUTINE trc_sms_fabm
498#endif
499
500   !!======================================================================
501END MODULE trcsms_fabm
Note: See TracBrowser for help on using the repository browser.