source: branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 9813

Last change on this file since 9813 was 9813, checked in by antsia, 2 years ago

delete iscplhsb, add iscpldiv and make the code readable

File size: 17.7 KB
Line 
1MODULE trasbc
2   !!==============================================================================
3   !!                       ***  MODULE  trasbc  ***
4   !! Ocean active tracers:  surface boundary condition
5   !!==============================================================================
6   !! History :  OPA  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
7   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
8   !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
10   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_sbc      : update the tracer trend at ocean surface
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and active tracers
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE dom_oce         ! ocean space domain variables
20   USE phycst          ! physical constant
21   USE sbcmod          ! ln_rnf 
22   USE sbcrnf          ! River runoff 
23   USE sbcisf          ! Ice shelf   
24   USE traqsr          ! solar radiation penetration
25   USE trd_oce         ! trends: ocean variables
26   USE trdtra          ! trends manager: tracers
27   !
28   USE in_out_manager  ! I/O manager
29   USE prtctl          ! Print control
30   USE iom
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory Allocation
33   USE timing          ! Timing
34   USE eosbn2
35   USE iscpldiv
36#if defined key_asminc   
37   USE asminc          ! Assimilation increment
38#endif
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   tra_sbc    ! routine called by step.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE tra_sbc ( kt )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE tra_sbc  ***
58      !!                   
59      !! ** Purpose :   Compute the tracer surface boundary condition trend of
60      !!      (flux through the interface, concentration/dilution effect)
61      !!      and add it to the general trend of tracer equations.
62      !!
63      !! ** Method :
64      !!      Following Roullet and Madec (2000), the air-sea flux can be divided
65      !!      into three effects: (1) Fext, external forcing;
66      !!      (2) Fwi, concentration/dilution effect due to water exchanged
67      !!         at the surface by evaporation, precipitations and runoff (E-P-R);
68      !!      (3) Fwe, tracer carried with the water that is exchanged.
69      !!            - salinity    : salt flux only due to freezing/melting
70      !!            sa = sa +  sfx / rau0 / e3t  for k=1
71      !!
72      !!      Fext, flux through the air-sea interface for temperature and salt:
73      !!            - temperature : heat flux q (w/m2). If penetrative solar
74      !!         radiation q is only the non solar part of the heat flux, the
75      !!         solar part is added in traqsr.F routine.
76      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
77      !!            - salinity    : no salt flux
78      !!
79      !!      The formulation for Fwb and Fwi vary according to the free
80      !!      surface formulation (linear or variable volume).
81      !!      * Linear free surface
82      !!            The surface freshwater flux modifies the ocean volume
83      !!         and thus the concentration of a tracer and the temperature.
84      !!         First order of the effect of surface freshwater exchange
85      !!         for salinity, it can be neglected on temperature (especially
86      !!         as the temperature of precipitations and runoffs is usually
87      !!         unknown).
88      !!            - temperature : we assume that the temperature of both
89      !!         precipitations and runoffs is equal to the SST, thus there
90      !!         is no additional flux since in this case, the concentration
91      !!         dilution effect is balanced by the net heat flux associated
92      !!         to the freshwater exchange (Fwe+Fwi=0):
93      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
94      !!            - salinity    : evaporation, precipitation and runoff
95      !!         water has a zero salinity  but there is a salt flux due to
96      !!         freezing/melting, thus:
97      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
98      !!                    + sfx    / rau0 / e3t
99      !!         where emp, the surface freshwater budget (evaporation minus
100      !!         precipitation minus runoff) given in kg/m2/s is divided
101      !!         by rau0 (density of sea water) to obtain m/s.   
102      !!         Note: even though Fwe does not appear explicitly for
103      !!         temperature in this routine, the heat carried by the water
104      !!         exchanged through the surface is part of the total heat flux
105      !!         forcing and must be taken into account in the global heat
106      !!         balance).
107      !!      * nonlinear free surface (variable volume, lk_vvl)
108      !!         contrary to the linear free surface case, Fwi is properly
109      !!         taken into account by using the true layer thicknesses to       
110      !!         calculate tracer content and advection. There is no need to
111      !!         deal with it in this routine.
112      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
113      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
114      !!
115      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
116      !!                with the tracer surface boundary condition
117      !!              - send trends to trdtra module (l_trdtra=T)
118      !!----------------------------------------------------------------------
119      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
120      !!
121      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
122      INTEGER  ::   ikt, ikb 
123      INTEGER  ::   nk_isf
124      REAL(wp) ::   zfact, z1_e3t, zdep
125      REAL(wp) ::   zalpha, zhk
126      REAL(wp) ::  zt_frz, zpress
127      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
128      !!----------------------------------------------------------------------
129      !
130      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
131      !
132      IF( kt == nit000 ) THEN
133         IF(lwp) WRITE(numout,*)
134         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
135         IF(lwp) WRITE(numout,*) '~~~~~~~ '
136      ENDIF
137
138      IF( l_trdtra ) THEN                    !* Save ta and sa trends
139         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
140         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
141         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
142      ENDIF
143
144!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
145      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
146         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
147         qsr(:,:) = 0.e0                     ! qsr set to zero
148      ENDIF
149
150      !----------------------------------------
151      !        EMP, SFX and QNS effects
152      !----------------------------------------
153      !                                          Set before sbc tracer content fields
154      !                                          ************************************
155      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
156         !                                         ! -----------------------------------
157         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
158              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
159            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
160            zfact = 0.5_wp
161            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
162            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
163         ELSE                                         ! No restart or restart not found: Euler forward time stepping
164            zfact = 1._wp
165            sbc_tsc(:,:,:) = 0._wp
166            sbc_tsc_b(:,:,:) = 0._wp
167         ENDIF
168      ELSE                                         ! Swap of forcing fields
169         !                                         ! ----------------------
170         zfact = 0.5_wp
171         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
172      ENDIF
173      !                                          Compute now sbc tracer content fields
174      !                                          *************************************
175
176                                                   ! Concentration dilution effect on (t,s) due to 
177                                                   ! evaporation, precipitation and qns, but not river runoff
178                                               
179      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
180         DO jj = 1, jpj
181            DO ji = 1, jpi 
182               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
183               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
184            END DO
185         END DO
186      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
187         DO jj = 2, jpj
188            DO ji = fs_2, fs_jpim1   ! vector opt.
189               ! temperature : heat flux
190               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
191                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
192               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
193               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
194                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
195            END DO
196         END DO
197         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
198         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
199      ENDIF
200      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff 
201      DO jn = 1, jpts
202         DO jj = 2, jpj
203            DO ji = fs_2, fs_jpim1   ! vector opt.
204               z1_e3t = zfact / fse3t(ji,jj,1)
205               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
206            END DO
207         END DO
208      END DO
209      !                                          Write in the ocean restart file
210      !                                          *******************************
211      IF( lrst_oce ) THEN
212         IF(lwp) WRITE(numout,*)
213         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
214            &                    'at it= ', kt,' date= ', ndastp
215         IF(lwp) WRITE(numout,*) '~~~~'
216         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
217         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
218         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
219         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
220      ENDIF
221      !
222      !
223
224      !---start the tendency with tracer source/sink from divergence correction
225      IF ( ln_iscpl .AND. kt == nit000) THEN       
226         CALL iscpl_div_tra( tsa )
227      ENDIF
228
229      !----------------------------------------
230      !       Ice Shelf effects (ISF)
231      !     tbl treated as in Losh (2008) JGR
232      !----------------------------------------
233      !
234      IF( nn_isf > 0 ) THEN
235         zfact = 0.5e0
236         DO jj = 2, jpj
237            DO ji = fs_2, fs_jpim1
238         
239               ikt = misfkt(ji,jj)
240               ikb = misfkb(ji,jj)
241   
242               ! level fully include in the ice shelf boundary layer
243               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
244               ! sign - because fwf sign of evapo (rnf sign of precip)
245               DO jk = ikt, ikb - 1
246               ! compute tfreez for the temperature correction (we add water at freezing temperature)
247               ! compute trend
248                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
249                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj)
250                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
251                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
252               END DO
253   
254               ! level partially include in ice shelf boundary layer
255               ! compute tfreez for the temperature correction (we add water at freezing temperature)
256               ! compute trend
257               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
258                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
259               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
260                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
261            END DO
262         END DO
263         IF( lrst_oce ) THEN
264            IF(lwp) WRITE(numout,*)
265            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
266               &                    'at it= ', kt,' date= ', ndastp
267            IF(lwp) WRITE(numout,*) '~~~~'
268            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
269            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
270            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
271            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
272            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
273         ENDIF
274      END IF
275      !
276      !----------------------------------------
277      !        River Runoff effects
278      !----------------------------------------
279      !
280      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
281         zfact = 0.5_wp
282         DO jj = 2, jpj 
283            DO ji = fs_2, fs_jpim1
284               IF( rnf(ji,jj) /= 0._wp ) THEN
285                  zdep = zfact / h_rnf(ji,jj)
286                  DO jk = 1, nk_rnf(ji,jj)
287                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
288                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
289                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
290                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
291                  END DO
292               ENDIF
293            END DO 
294         END DO 
295      ENDIF
296
297      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
298      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
299
300#if defined key_asminc
301! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
302! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
303! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
304      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
305         DO jj = 2, jpj 
306            DO ji = fs_2, fs_jpim1
307               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
308               DO jk = 1, jpkm1
309                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
310                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
311                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
312                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
313               END DO
314            END DO 
315         END DO 
316      ENDIF
317#endif
318 
319      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
320         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
321         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
322         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
323         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
324         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
325      ENDIF
326      !
327      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
328         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
329      !
330      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
331      !
332   END SUBROUTINE tra_sbc
333
334   !!======================================================================
335END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.