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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 13.2 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, Kmm, Krhs )
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      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
76      !
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
80      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
81      !!----------------------------------------------------------------------
82      !
83      IF( ln_timing )   CALL timing_start('tra_sbc')
84      !
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
90      !
91      IF( l_trdtra ) THEN                    !* Save ta and sa trends
92         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
93         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
94         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
95      ENDIF
96      !
97!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
98      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
99         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
100         qsr(:,:) = 0._wp                     ! qsr set to zero
101      ENDIF
102
103      !----------------------------------------
104      !        EMP, SFX and QNS effects
105      !----------------------------------------
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
109              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
110            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
111            zfact = 0.5_wp
112            sbc_tsc(:,:,:) = 0._wp
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
115         ELSE                                   ! No restart or restart not found: Euler forward time stepping
116            zfact = 1._wp
117            sbc_tsc(:,:,:) = 0._wp
118            sbc_tsc_b(:,:,:) = 0._wp
119         ENDIF
120      ELSE                                !* other time-steps: swap of forcing fields
121         zfact = 0.5_wp
122         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
123      ENDIF
124      !                             !==  Now sbc tracer content fields  ==!
125      DO jj = 2, jpj
126         DO ji = fs_2, fs_jpim1   ! vector opt.
127            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
128            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
129         END DO
130      END DO
131      IF( ln_linssh ) THEN                !* linear free surface 
132         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
133            DO ji = fs_2, fs_jpim1   ! vector opt.
134               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
135               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)
136            END DO
137         END DO                                 !==>> output c./d. term
138         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )
139         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
140      ENDIF
141      !
142      DO jn = 1, jpts               !==  update tracer trend  ==!
143         DO jj = 2, jpj
144            DO ji = fs_2, fs_jpim1   ! vector opt. 
145               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)
146            END DO
147         END DO
148      END DO
149      !                 
150      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
151         IF( lwxios ) CALL iom_swap(      cwxios_context          )
152         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem), ldxios = lwxios )
153         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal), ldxios = lwxios )
154         IF( lwxios ) CALL iom_swap(      cxios_context          )
155      ENDIF
156      !
157      !----------------------------------------
158      !       Ice Shelf effects (ISF)
159      !     tbl treated as in Losh (2008) JGR
160      !----------------------------------------
161      !
162!!gm BUG ?   Why no differences between non-linear and linear free surface ?
163!!gm         probably taken into account in r1_hisf_tbl : to be verified
164      IF( ln_isf ) THEN
165         zfact = 0.5_wp
166         DO jj = 2, jpj
167            DO ji = fs_2, fs_jpim1
168               !
169               ikt = misfkt(ji,jj)
170               ikb = misfkb(ji,jj)
171               !
172               ! level fully include in the ice shelf boundary layer
173               ! sign - because fwf sign of evapo (rnf sign of precip)
174               DO jk = ikt, ikb - 1
175               ! compute trend
176                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
177                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
178                     &           * r1_hisf_tbl(ji,jj)
179               END DO
180   
181               ! level partially include in ice shelf boundary layer
182               ! compute trend
183               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
184                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
185                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
186
187            END DO
188         END DO
189      END IF
190      !
191      !----------------------------------------
192      !        River Runoff effects
193      !----------------------------------------
194      !
195      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
196         zfact = 0.5_wp
197         DO jj = 2, jpj 
198            DO ji = fs_2, fs_jpim1
199               IF( rnf(ji,jj) /= 0._wp ) THEN
200                  zdep = zfact / h_rnf(ji,jj)
201                  DO jk = 1, nk_rnf(ji,jj)
202                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 &
203                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
204                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 &
205                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
206                  END DO
207               ENDIF
208            END DO 
209         END DO 
210      ENDIF
211
212      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
213      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
214
215#if defined key_asminc
216      !
217      !----------------------------------------
218      !        Assmilation effects
219      !----------------------------------------
220      !
221      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
222          !
223         IF( ln_linssh ) THEN
224            DO jj = 2, jpj 
225               DO ji = fs_2, fs_jpim1
226                  ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1)
227                  tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim
228                  tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim
229               END DO
230            END DO
231         ELSE
232            DO jj = 2, jpj 
233               DO ji = fs_2, fs_jpim1
234                  ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
235                  tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim
236                  tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim
237               END DO 
238            END DO 
239         ENDIF
240         !
241      ENDIF
242      !
243#endif
244      !
245      !----------------------------------------
246      !        Ice Sheet coupling imbalance correction to have conservation
247      !----------------------------------------
248      !
249      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff
250         DO jk = 1,jpk
251            DO jj = 2, jpj 
252               DO ji = fs_2, fs_jpim1
253                  zdep = 1._wp / e3t_n(ji,jj,jk) 
254                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
255                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep 
256               END DO 
257            END DO 
258         END DO
259      ENDIF
260
261      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
262         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
263         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
264         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
265         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
266         DEALLOCATE( ztrdt , ztrds ) 
267      ENDIF
268      !
269      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
270         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
271      !
272      IF( ln_timing )   CALL timing_stop('tra_sbc')
273      !
274   END SUBROUTINE tra_sbc
275
276   !!======================================================================
277END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.