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.
trasbc.F90 in NEMO/branches/UKMO/dev_10448_WAD_SBC_BUGFIX/src/OCE/TRA – NEMO

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, 5 years 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
RevLine 
[3]1MODULE trasbc
2   !!==============================================================================
3   !!                       ***  MODULE  trasbc  ***
4   !! Ocean active tracers:  surface boundary condition
5   !!==============================================================================
[2528]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
[5120]11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing
[3]12   !!----------------------------------------------------------------------
[503]13
14   !!----------------------------------------------------------------------
[6140]15   !!   tra_sbc       : update the tracer trend at ocean surface
[3]16   !!----------------------------------------------------------------------
[6140]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
[9023]29#if defined key_asminc   
30   USE asminc         ! Assimilation increment
31#endif
[4990]32   !
[6140]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
[3]38
39   IMPLICIT NONE
40   PRIVATE
41
[6140]42   PUBLIC   tra_sbc   ! routine called by step.F90
[3]43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
[9598]47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]48   !! $Id$
[10068]49   !! Software governed by the CeCILL license (see ./LICENSE)
[3]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      !!
[6140]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.
[664]70      !!
[6140]71      !! ** Action  : - Update tsa with the surface boundary condition trend
72      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
[503]73      !!----------------------------------------------------------------------
[2528]74      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[6140]75      !
[9023]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
[9019]79      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
[3]80      !!----------------------------------------------------------------------
[3294]81      !
[9019]82      IF( ln_timing )   CALL timing_start('tra_sbc')
[3294]83      !
[3]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
[6140]89      !
[4990]90      IF( l_trdtra ) THEN                    !* Save ta and sa trends
[9019]91         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
[7753]92         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
93         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]94      ENDIF
[6140]95      !
96!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
[2528]97      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
[7753]98         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
99         qsr(:,:) = 0._wp                     ! qsr set to zero
[2528]100      ENDIF
[3]101
[2528]102      !----------------------------------------
[4990]103      !        EMP, SFX and QNS effects
[2528]104      !----------------------------------------
[6140]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
[2528]108              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
[6140]109            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
[4990]110            zfact = 0.5_wp
[7753]111            sbc_tsc(:,:,:) = 0._wp
[9367]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
[6140]114         ELSE                                   ! No restart or restart not found: Euler forward time stepping
[4990]115            zfact = 1._wp
[7753]116            sbc_tsc(:,:,:) = 0._wp
117            sbc_tsc_b(:,:,:) = 0._wp
[2528]118         ENDIF
[6140]119      ELSE                                !* other time-steps: swap of forcing fields
[4990]120         zfact = 0.5_wp
[7753]121         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
[2528]122      ENDIF
[6140]123      !                             !==  Now sbc tracer content fields  ==!
124      DO jj = 2, jpj
125         DO ji = fs_2, fs_jpim1   ! vector opt.
[10456]126            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
[6140]127            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
[3]128         END DO
[6140]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
[2528]132            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]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)
[2528]135            END DO
[6140]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) )
[2528]139      ENDIF
[6140]140      !
141      DO jn = 1, jpts               !==  update tracer trend  ==!
[2528]142         DO jj = 2, jpj
[6140]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)
[2528]145            END DO
146         END DO
[3]147      END DO
[6140]148      !                 
149      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
[9367]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          )
[2528]154      ENDIF
155      !
156      !----------------------------------------
[4990]157      !       Ice Shelf effects (ISF)
158      !     tbl treated as in Losh (2008) JGR
159      !----------------------------------------
160      !
[6140]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
[4990]165         DO jj = 2, jpj
166            DO ji = fs_2, fs_jpim1
[6140]167               !
[4990]168               ikt = misfkt(ji,jj)
169               ikb = misfkb(ji,jj)
[6140]170               !
[4990]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
[6140]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)
[4990]178               END DO
179   
180               ! level partially include in ice shelf boundary layer
181               ! compute trend
[6140]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
[4990]186            END DO
187         END DO
188      END IF
189      !
190      !----------------------------------------
[2528]191      !        River Runoff effects
192      !----------------------------------------
193      !
[3764]194      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
195         zfact = 0.5_wp
[2528]196         DO jj = 2, jpj 
197            DO ji = fs_2, fs_jpim1
[3764]198               IF( rnf(ji,jj) /= 0._wp ) THEN
199                  zdep = zfact / h_rnf(ji,jj)
[2528]200                  DO jk = 1, nk_rnf(ji,jj)
[6140]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 
[2715]205                  END DO
[2528]206               ENDIF
[2715]207            END DO 
208         END DO 
[3764]209      ENDIF
[6472]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
[9023]214#if defined key_asminc
[6140]215      !
216      !----------------------------------------
[9023]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      !----------------------------------------
[6140]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) 
[9124]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 
[6140]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
[7753]261         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
262         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[4990]263         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
264         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
[9019]265         DEALLOCATE( ztrdt , ztrds ) 
[216]266      ENDIF
[503]267      !
[2528]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' )
[503]270      !
[9019]271      IF( ln_timing )   CALL timing_stop('tra_sbc')
[3294]272      !
[3]273   END SUBROUTINE tra_sbc
274
275   !!======================================================================
276END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.