Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
MODULE trasbc
!!==============================================================================
!! *** MODULE trasbc ***
!! Ocean active tracers: surface boundary condition
!!==============================================================================
!! History : OPA ! 1998-10 (G. Madec, G. Roullet, M. Imbard) Original code
!! 8.2 ! 2001-02 (D. Ludicone) sea ice and free surface
!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
!! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps
!! - ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC
!! 3.6 ! 2014-11 (P. Mathiot) isf melting forcing
!! 4.1 ! 2019-09 (P. Mathiot) isf moved in traisf
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! tra_sbc : update the tracer trend at ocean surface
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and active tracers
USE sbc_oce ! surface boundary condition: ocean
USE dom_oce ! ocean space domain variables
USE phycst ! physical constant
USE eosbn2 ! Equation Of State
USE sbcmod ! ln_rnf
USE sbcrnf ! River runoff
USE traqsr ! solar radiation penetration
USE trd_oce ! trends: ocean variables
USE trdtra ! trends manager: tracers
#if defined key_asminc
USE asminc ! Assimilation increment
#endif
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
USE iom ! xIOS server
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC tra_sbc ! routine called by step.F90
PUBLIC tra_sbc_RK3 ! routine called by stprk3_stg.F90
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trasbc.F90 15379 2021-10-15 09:05:45Z techene $
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE tra_sbc ( kt, Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_sbc ***
!!
!! ** Purpose : Compute the tracer surface boundary condition trend of
!! (flux through the interface, concentration/dilution effect)
!! and add it to the general trend of tracer equations.
!!
!! ** Method : The (air+ice)-sea flux has two components:
!! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
!! (2) Fwe , tracer carried with the water that is exchanged with air+ice.
!! The input forcing fields (emp, rnf, sfx) contain Fext+Fwe,
!! they are simply added to the tracer trend (ts(Krhs)).
!! In linear free surface case (ln_linssh=T), the volume of the
!! ocean does not change with the water exchanges at the (air+ice)-sea
!! interface. Therefore another term has to be added, to mimic the
!! concentration/dilution effect associated with water exchanges.
!!
!! ** Action : - Update ts(Krhs) with the surface boundary condition trend
!! - send trends to trdtra module for further diagnostics(l_trdtra=T)
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt ! ocean time-step index
INTEGER, INTENT(in ) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer Eq.
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ikt, ikb ! local integers
REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_sbc')
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
ENDIF
ENDIF
!
IF( l_trdtra ) THEN !* Save ta and sa trends
ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
ENDIF
!
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
ENDIF
!----------------------------------------
! EMP, SFX and QNS effects
!----------------------------------------
! !== Set before sbc tracer content fields ==!
IF( kt == nit000 ) THEN !* 1st time-step
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN ! Restart: read in restart file
zfact = 0.5_wp
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF(lwp) WRITE(numout,*) ' nit000-1 sbc tracer content field read in the restart file'
sbc_tsc(:,:,:) = 0._wp
CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend
CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend
ENDIF
ELSE ! No restart or restart not found: Euler forward time stepping
zfact = 1._wp
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
sbc_tsc(ji,jj,:) = 0._wp
sbc_tsc_b(ji,jj,:) = 0._wp
END_2D
ENDIF
ELSE !* other time-steps: swap of forcing fields
zfact = 0.5_wp
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
END_2D
ENDIF
! !== Now sbc tracer content fields ==!
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux
sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting
END_2D
IF( ln_linssh ) THEN !* linear free surface
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) !==>> add concentration/dilution effect due to constant volume cell
sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
END_2D !==>> output c./d. term
IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
ENDIF
!
DO jn = 1, jpts !== update tracer trend ==!
DO_2D( 0, 0, 0, 0 )
pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) &
& / e3t(ji,jj,1,Kmm)
END_2D
END DO
!
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==!
CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
ENDIF
ENDIF
!
!----------------------------------------
! River Runoff effects
!----------------------------------------
!
IF( ln_rnf ) THEN ! input of heat and salt due to river runoff
zfact = 0.5_wp
DO_2D( 0, 0, 0, 0 )
IF( rnf(ji,jj) /= 0._wp ) THEN
zdep = zfact / h_rnf(ji,jj)
DO jk = 1, nk_rnf(ji,jj)
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
& + ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
IF( ln_rnf_sal ) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) &
& + ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep
END DO
ENDIF
END_2D
ENDIF
IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst
IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss
#if defined key_asminc
!
!----------------------------------------
! Assmilation effects
!----------------------------------------
!
IF( ln_sshinc ) THEN ! input of heat and salt due to assimilation
!
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
END_2D
ENDIF
!
ENDIF
!
#endif
!
IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
DEALLOCATE( ztrdt , ztrds )
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc - Ta: ', mask1=tmask, &
& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
IF( ln_timing ) CALL timing_stop('tra_sbc')
!
END SUBROUTINE tra_sbc
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
SUBROUTINE tra_sbc_RK3 ( kt, Kmm, pts, Krhs, kstg )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_sbc_RK3 ***
!!
!! ** Purpose : Compute the tracer surface boundary condition trend of
!! (flux through the interface, concentration/dilution effect)
!! and add it to the general trend of tracer equations.
!!
!! ** Method : The (air+ice)-sea flux has two components:
!! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
!! (2) Fwe , tracer carried with the water that is exchanged with air+ice.
!! The input forcing fields (emp, rnf, sfx) contain Fext+Fwe,
!! they are simply added to the tracer trend (ts(Krhs)).
!! In linear free surface case (ln_linssh=T), the volume of the
!! ocean does not change with the water exchanges at the (air+ice)-sea
!! interface. Therefore another term has to be added, to mimic the
!! concentration/dilution effect associated with water exchanges.
!!
!! ** Action : - Update ts(Krhs) with the surface boundary condition trend
!! - send trends to trdtra module for further diagnostics(l_trdtra=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices
INTEGER , INTENT(in ) :: kstg ! RK3 stage index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer Eq.
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: z1_rho0_e3t, zdep, ztim ! local scalar
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_sbc_RK3')
!
IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_sbc_RK3 : TRAcer Surface Boundary Condition'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
ENDIF
!
IF( l_trdtra ) THEN !* Save ta and sa trends
ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
ENDIF
!
ENDIF
!
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr .AND. kstg == 1) THEN ! no solar radiation penetration
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
ENDIF
!----------------------------------------
! EMP, SFX and QNS effects
!----------------------------------------
! !== update tracer trend ==!
SELECT CASE( kstg )
!
CASE( 1 , 2 ) != stage 1 and 2 =! only in non linear ssh
!
IF( .NOT.ln_linssh ) THEN !* only heat and salt fluxes associated with mass fluxes
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm) * z1_rho0_e3t
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm) * z1_rho0_e3t
END_2D
ENDIF
!
CASE( 3 )
!
IF( ln_linssh ) THEN !* linear free surface
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + ( r1_rcp * qns(ji,jj) & ! non solar heat flux
& + emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm) ) * z1_rho0_e3t ! add concentration/dilution effect due to constant volume cell
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + ( sfx(ji,jj) & ! salt flux due to freezing/melting
& + emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm) ) * z1_rho0_e3t ! add concentration/dilution effect due to constant volume cell
END_2D
IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
ENDIF
ELSE
DO_2D( 0, 0, 0, 0 )
z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + r1_rcp * qns(ji,jj) * z1_rho0_e3t
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + sfx(ji,jj) * z1_rho0_e3t
END_2D
ENDIF
END SELECT
!
!
!----------------------------------------
! River Runoff effects
!----------------------------------------
!
IF( ln_rnf ) THEN ! input of heat and salt due to river runoff
DO_2D( 0, 0, 0, 0 )
IF( rnf(ji,jj) /= 0._wp ) THEN
zdep = 1._wp / h_rnf(ji,jj)
DO jk = 1, nk_rnf(ji,jj)
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + rnf_tsc(ji,jj,jp_tem) * zdep
IF( ln_rnf_sal ) pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + rnf_tsc(ji,jj,jp_sal) * zdep
END DO
ENDIF
END_2D
ENDIF
!
IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst
IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss
ENDIF
#if defined key_asminc
!
!----------------------------------------
! Assmilation effects
!----------------------------------------
!
IF( ln_sshinc .AND. kstg == 3 ) THEN ! input of heat and salt due to assimilation
!
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
END_2D
ENDIF
!
ENDIF
!
#endif
!
IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics
IF( ntile == 0 .OR. ntile == nijtile ) THEN
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
DEALLOCATE( ztrdt , ztrds )
ENDIF
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc - Ta: ', mask1=tmask, &
& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
!
IF( ln_timing ) CALL timing_stop('tra_sbc_RK3')
!
END SUBROUTINE tra_sbc_RK3
!!======================================================================
END MODULE trasbc