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 @ 14072

Last change on this file since 14072 was 14072, checked in by laurent, 3 years ago

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
  • Property svn:keywords set to Id
File size: 12.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   !!            4.1  !  2019-09  (P. Mathiot) isf moved in traisf
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_sbc       : update the tracer trend at ocean surface
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE sbc_oce        ! surface boundary condition: ocean
20   USE dom_oce        ! ocean space domain variables
21   USE phycst         ! physical constant
22   USE eosbn2         ! Equation Of State
23   USE sbcmod         ! ln_rnf
24   USE sbcrnf         ! River runoff
25   USE traqsr         ! solar radiation penetration
26   USE trd_oce        ! trends: ocean variables
27   USE trdtra         ! trends manager: tracers
28#if defined key_asminc
29   USE asminc         ! Assimilation increment
30#endif
31   !
32   USE in_out_manager ! I/O manager
33   USE prtctl         ! Print control
34   USE iom            ! xIOS server
35   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
36   USE timing         ! Timing
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   tra_sbc   ! routine called by step.F90
42
43   !! * Substitutions
44#  include "do_loop_substitute.h90"
45#  include "domzgr_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, pts, 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) contain Fext+Fwe,
65      !!             they are simply added to the tracer trend (ts(Krhs)).
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 ts(Krhs) 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      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts        ! active tracers and RHS of tracer Eq.
77      !
78      INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices
79      INTEGER  ::   ikt, ikb, isi, iei, isj, iej ! local integers
80      REAL(wp) ::   zfact, z1_e3t, zdep, ztim    ! local scalar
81      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
82      !!----------------------------------------------------------------------
83      !
84      IF( ln_timing )   CALL timing_start('tra_sbc')
85      !
86      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
87         IF( kt == nit000 ) THEN
88            IF(lwp) WRITE(numout,*)
89            IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
90            IF(lwp) WRITE(numout,*) '~~~~~~~ '
91         ENDIF
92      ENDIF
93      !
94      IF( l_trdtra ) THEN                    !* Save ta and sa trends
95         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
96         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
97         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
98      ENDIF
99      !
100      ! NOTE: [tiling-comms-merge] Many DO loop bounds changed to avoid changing results when using tiling. Some bounds were also adjusted to account for those changed in tra_atf
101      IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
102      IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF
103      IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF
104      IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF
105
106!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
107      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
108         DO_2D( isj, iej, isi, iei )
109            qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns
110            qsr(ji,jj) = 0._wp                        ! qsr set to zero
111         END_2D
112      ENDIF
113
114      !----------------------------------------
115      !        EMP, SFX and QNS effects
116      !----------------------------------------
117      !                             !==  Set before sbc tracer content fields  ==!
118      IF( kt == nit000 ) THEN             !* 1st time-step
119         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN      ! Restart: read in restart file
120            zfact = 0.5_wp
121            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
122               IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
123               sbc_tsc(:,:,:) = 0._wp
124               CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
125               CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
126            ENDIF
127         ELSE                                             ! No restart or restart not found: Euler forward time stepping
128            zfact = 1._wp
129            DO_2D( isj, iej, isi, iei )
130               sbc_tsc(ji,jj,:) = 0._wp
131               sbc_tsc_b(ji,jj,:) = 0._wp
132            END_2D
133         ENDIF
134      ELSE                                !* other time-steps: swap of forcing fields
135         zfact = 0.5_wp
136         DO_2D( isj, iej, isi, iei )
137            sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
138         END_2D
139      ENDIF
140      !                             !==  Now sbc tracer content fields  ==!
141      DO_2D( isj, iej, isi, iei )
142         sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
143         sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
144      END_2D
145      IF( ln_linssh ) THEN                !* linear free surface
146         DO_2D( isj, iej, isi, iei )                    !==>> add concentration/dilution effect due to constant volume cell
147            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
148            sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
149         END_2D                                 !==>> output c./d. term
150         IF( ntile == 0 .OR. ntile == nijtile )  THEN             ! Do only on the last tile
151            IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
152            IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
153         ENDIF
154      ENDIF
155      !
156      DO jn = 1, jpts               !==  update tracer trend  ==!
157         ! NOTE: [tiling-comms-merge] This looped over nn_hls, which changes the results when using tiling
158         DO_2D( 0, 0, 0, 0 )
159            pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) )    &
160               &                                                / e3t(ji,jj,1,Kmm)
161         END_2D
162      END DO
163      !
164      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
165         IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
166            CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
167            CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
168         ENDIF
169      ENDIF
170      !
171      !----------------------------------------
172      !        River Runoff effects
173      !----------------------------------------
174      !
175      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
176         zfact = 0.5_wp
177         DO_2D( 0, 0, 0, 0 )
178            IF( rnf(ji,jj) /= 0._wp ) THEN
179               zdep = zfact / h_rnf(ji,jj)
180               DO jk = 1, nk_rnf(ji,jj)
181                                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)                                  &
182                                        &                      +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
183                  IF( ln_rnf_sal )   pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)                                  &
184                                        &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep
185               END DO
186            ENDIF
187         END_2D
188      ENDIF
189
190      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
191         IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
192         IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
193      ENDIF
194
195#if defined key_asminc
196      !
197      !----------------------------------------
198      !        Assmilation effects
199      !----------------------------------------
200      !
201      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
202          !
203         IF( ln_linssh ) THEN
204            DO_2D( 0, 0, 0, 0 )
205               ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
206               pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
207               pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim
208            END_2D
209         ELSE
210            DO_2D( 0, 0, 0, 0 )
211               ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
212               pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
213               pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
214            END_2D
215         ENDIF
216         !
217      ENDIF
218      !
219#endif
220      !
221      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
222         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
223         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
224         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
225         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
226         DEALLOCATE( ztrdt , ztrds )
227      ENDIF
228      !
229      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
230         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
231      !
232      IF( ln_timing )   CALL timing_stop('tra_sbc')
233      !
234   END SUBROUTINE tra_sbc
235
236   !!======================================================================
237END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.