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

source: branches/UKMO/dev_r5518_nemo_fabm_ukmo/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90 @ 7829

Last change on this file since 7829 was 7829, checked in by dford, 7 years ago

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