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/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/trasbc.F90 @ 10068

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

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