source: NEMO/branches/UKMO/dev_10448_WAD_SBC_BUGFIX/src/OCE/TRA/trasbc.F90 @ 10456

Last change on this file since 10456 was 10456, checked in by deazer, 22 months ago

Added option to taper sbc fluxes near very shallow water when using WAD
Corrected some IO bugs in dia25h, diatmb for WAD case.
User has control of the tapering. At what depth to start it, and at what fraction to start
the tanh tapering. At the WAD limit SBC is turned off completely.
Dry cells do not have any communication with the atmosphere
To DO: Documentation update.
Although not all sette tests are passed (AGRIF etc.)
it does no worse than the trunk at the revision the branch is made

  • Property svn:keywords set to Id
File size: 13.1 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 eosbn2         ! Equation Of State
22   USE sbcmod         ! ln_rnf 
23   USE sbcrnf         ! River runoff 
24   USE sbcisf         ! Ice shelf   
25   USE iscplini       ! Ice sheet coupling
26   USE traqsr         ! solar radiation penetration
27   USE trd_oce        ! trends: ocean variables
28   USE trdtra         ! trends manager: tracers
29#if defined key_asminc   
30   USE asminc         ! Assimilation increment
31#endif
32   !
33   USE in_out_manager ! I/O manager
34   USE prtctl         ! Print control
35   USE iom            ! xIOS server
36   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
37   USE timing         ! Timing
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   tra_sbc   ! routine called by step.F90
43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE tra_sbc ( kt )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_sbc  ***
56      !!                   
57      !! ** Purpose :   Compute the tracer surface boundary condition trend of
58      !!      (flux through the interface, concentration/dilution effect)
59      !!      and add it to the general trend of tracer equations.
60      !!
61      !! ** Method :   The (air+ice)-sea flux has two components:
62      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
63      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.
64      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe,
65      !!             they are simply added to the tracer trend (tsa).
66      !!               In linear free surface case (ln_linssh=T), the volume of the
67      !!             ocean does not change with the water exchanges at the (air+ice)-sea
68      !!             interface. Therefore another term has to be added, to mimic the
69      !!             concentration/dilution effect associated with water exchanges.
70      !!
71      !! ** Action  : - Update tsa with the surface boundary condition trend
72      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
73      !!----------------------------------------------------------------------
74      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
75      !
76      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
77      INTEGER  ::   ikt, ikb                    ! local integers
78      REAL(wp) ::   zfact, z1_e3t, zdep, ztim   ! local scalar
79      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
80      !!----------------------------------------------------------------------
81      !
82      IF( ln_timing )   CALL timing_start('tra_sbc')
83      !
84      IF( kt == nit000 ) THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
87         IF(lwp) WRITE(numout,*) '~~~~~~~ '
88      ENDIF
89      !
90      IF( l_trdtra ) THEN                    !* Save ta and sa trends
91         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
92         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
93         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
94      ENDIF
95      !
96!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
97      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
98         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
99         qsr(:,:) = 0._wp                     ! qsr set to zero
100      ENDIF
101
102      !----------------------------------------
103      !        EMP, SFX and QNS effects
104      !----------------------------------------
105      !                             !==  Set before sbc tracer content fields  ==!
106      IF( kt == nit000 ) THEN             !* 1st time-step
107         IF( ln_rstart .AND.    &               ! Restart: read in restart file
108              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
109            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
110            zfact = 0.5_wp
111            sbc_tsc(:,:,:) = 0._wp
112            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content sbc trend
113            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content sbc trend
114         ELSE                                   ! No restart or restart not found: Euler forward time stepping
115            zfact = 1._wp
116            sbc_tsc(:,:,:) = 0._wp
117            sbc_tsc_b(:,:,:) = 0._wp
118         ENDIF
119      ELSE                                !* other time-steps: swap of forcing fields
120         zfact = 0.5_wp
121         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
122      ENDIF
123      !                             !==  Now sbc tracer content fields  ==!
124      DO jj = 2, jpj
125         DO ji = fs_2, fs_jpim1   ! vector opt.
126            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
127            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
128         END DO
129      END DO
130      IF( ln_linssh ) THEN                !* linear free surface 
131         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
132            DO ji = fs_2, fs_jpim1   ! vector opt.
133               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
134               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)
135            END DO
136         END DO                                 !==>> output c./d. term
137         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )
138         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
139      ENDIF
140      !
141      DO jn = 1, jpts               !==  update tracer trend  ==!
142         DO jj = 2, jpj
143            DO ji = fs_2, fs_jpim1   ! vector opt. 
144               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1)
145            END DO
146         END DO
147      END DO
148      !                 
149      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
150         IF( lwxios ) CALL iom_swap(      cwxios_context          )
151         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem), ldxios = lwxios )
152         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal), ldxios = lwxios )
153         IF( lwxios ) CALL iom_swap(      cxios_context          )
154      ENDIF
155      !
156      !----------------------------------------
157      !       Ice Shelf effects (ISF)
158      !     tbl treated as in Losh (2008) JGR
159      !----------------------------------------
160      !
161!!gm BUG ?   Why no differences between non-linear and linear free surface ?
162!!gm         probably taken into account in r1_hisf_tbl : to be verified
163      IF( ln_isf ) THEN
164         zfact = 0.5_wp
165         DO jj = 2, jpj
166            DO ji = fs_2, fs_jpim1
167               !
168               ikt = misfkt(ji,jj)
169               ikb = misfkb(ji,jj)
170               !
171               ! level fully include in the ice shelf boundary layer
172               ! sign - because fwf sign of evapo (rnf sign of precip)
173               DO jk = ikt, ikb - 1
174               ! compute trend
175                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
176                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
177                     &           * r1_hisf_tbl(ji,jj)
178               END DO
179   
180               ! level partially include in ice shelf boundary layer
181               ! compute trend
182               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
183                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
184                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
185
186            END DO
187         END DO
188      END IF
189      !
190      !----------------------------------------
191      !        River Runoff effects
192      !----------------------------------------
193      !
194      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
195         zfact = 0.5_wp
196         DO jj = 2, jpj 
197            DO ji = fs_2, fs_jpim1
198               IF( rnf(ji,jj) /= 0._wp ) THEN
199                  zdep = zfact / h_rnf(ji,jj)
200                  DO jk = 1, nk_rnf(ji,jj)
201                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 &
202                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
203                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 &
204                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
205                  END DO
206               ENDIF
207            END DO 
208         END DO 
209      ENDIF
210
211      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
212      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
213
214#if defined key_asminc
215      !
216      !----------------------------------------
217      !        Assmilation effects
218      !----------------------------------------
219      !
220      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
221          !
222         IF( ln_linssh ) THEN
223            DO jj = 2, jpj 
224               DO ji = fs_2, fs_jpim1
225                  ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1)
226                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim
227                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim
228               END DO
229            END DO
230         ELSE
231            DO jj = 2, jpj 
232               DO ji = fs_2, fs_jpim1
233                  ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
234                  tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim
235                  tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim
236               END DO 
237            END DO 
238         ENDIF
239         !
240      ENDIF
241      !
242#endif
243      !
244      !----------------------------------------
245      !        Ice Sheet coupling imbalance correction to have conservation
246      !----------------------------------------
247      !
248      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff
249         DO jk = 1,jpk
250            DO jj = 2, jpj 
251               DO ji = fs_2, fs_jpim1
252                  zdep = 1._wp / e3t_n(ji,jj,jk) 
253                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
254                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep 
255               END DO 
256            END DO 
257         END DO
258      ENDIF
259
260      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
261         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
262         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
263         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
264         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
265         DEALLOCATE( ztrdt , ztrds ) 
266      ENDIF
267      !
268      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
269         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
270      !
271      IF( ln_timing )   CALL timing_stop('tra_sbc')
272      !
273   END SUBROUTINE tra_sbc
274
275   !!======================================================================
276END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.