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 branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 5866

Last change on this file since 5866 was 5866, checked in by gm, 8 years ago

#1613: vvl by default: add ln_linssh and remove key_vvl
  • Property svn:keywords set to Id
File size: 15.5 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 sbcmod          ! ln_rnf 
22   USE sbcrnf          ! River runoff 
23   USE sbcisf          ! Ice shelf   
24   USE traqsr          ! solar radiation penetration
25   USE trd_oce         ! trends: ocean variables
26   USE trdtra          ! trends manager: tracers
27   !
28   USE in_out_manager  ! I/O manager
29   USE prtctl          ! Print control
30   USE iom
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory Allocation
33   USE timing          ! Timing
34   USE eosbn2
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_sbc    ! routine called by step.F90
40
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE tra_sbc ( kt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_sbc  ***
53      !!                   
54      !! ** Purpose :   Compute the tracer surface boundary condition trend of
55      !!      (flux through the interface, concentration/dilution effect)
56      !!      and add it to the general trend of tracer equations.
57      !!
58      !! ** Method :
59      !!      Following Roullet and Madec (2000), the air-sea flux can be divided
60      !!      into three effects: (1) Fext, external forcing;
61      !!      (2) Fwi, concentration/dilution effect due to water exchanged
62      !!         at the surface by evaporation, precipitations and runoff (E-P-R);
63      !!      (3) Fwe, tracer carried with the water that is exchanged.
64      !!            - salinity    : salt flux only due to freezing/melting
65      !!            sa = sa +  sfx / rau0 / e3t  for k=1
66      !!
67      !!      Fext, flux through the air-sea interface for temperature and salt:
68      !!            - temperature : heat flux q (w/m2). If penetrative solar
69      !!         radiation q is only the non solar part of the heat flux, the
70      !!         solar part is added in traqsr.F routine.
71      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
72      !!            - salinity    : no salt flux
73      !!
74      !!      The formulation for Fwb and Fwi vary according to the free
75      !!      surface formulation (linear or variable volume).
76      !!      * Linear free surface (ln_linssh=T)
77      !!            The surface freshwater flux modifies the ocean volume
78      !!         and thus the concentration of a tracer and the temperature.
79      !!         First order of the effect of surface freshwater exchange
80      !!         for salinity, it can be neglected on temperature (especially
81      !!         as the temperature of precipitations and runoffs is usually
82      !!         unknown).
83      !!            - temperature : we assume that the temperature of both
84      !!         precipitations and runoffs is equal to the SST, thus there
85      !!         is no additional flux since in this case, the concentration
86      !!         dilution effect is balanced by the net heat flux associated
87      !!         to the freshwater exchange (Fwe+Fwi=0):
88      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
89      !!            - salinity    : evaporation, precipitation and runoff
90      !!         water has a zero salinity  but there is a salt flux due to
91      !!         freezing/melting, thus:
92      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
93      !!                    + sfx    / rau0 / e3t
94      !!         where emp, the surface freshwater budget (evaporation minus
95      !!         precipitation minus runoff) given in kg/m2/s is divided
96      !!         by rau0 (density of sea water) to obtain m/s.   
97      !!         Note: even though Fwe does not appear explicitly for
98      !!         temperature in this routine, the heat carried by the water
99      !!         exchanged through the surface is part of the total heat flux
100      !!         forcing and must be taken into account in the global heat
101      !!         balance).
102      !!      * nonlinear free surface (ln_linssh=F)
103      !!         contrary to the linear free surface case, Fwi is properly
104      !!         taken into account by using the true layer thicknesses to       
105      !!         calculate tracer content and advection. There is no need to
106      !!         deal with it in this routine.
107      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
108      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
109      !!
110      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
111      !!                with the tracer surface boundary condition
112      !!              - send trends to trdtra module (l_trdtra=T)
113      !!----------------------------------------------------------------------
114      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
115      !!
116      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
117      INTEGER  ::   ikt, ikb 
118      INTEGER  ::   nk_isf
119      REAL(wp) ::   zfact, z1_e3t, zdep
120      REAL(wp) ::   zalpha, zhk
121      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
122      !!----------------------------------------------------------------------
123      !
124      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
125      !
126      IF( kt == nit000 ) THEN
127         IF(lwp) WRITE(numout,*)
128         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
129         IF(lwp) WRITE(numout,*) '~~~~~~~ '
130      ENDIF
131
132      IF( l_trdtra ) THEN                    !* Save ta and sa trends
133         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
134         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
135         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
136      ENDIF
137
138!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
139      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
140         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
141         qsr(:,:) = 0.e0                     ! qsr set to zero
142      ENDIF
143
144      !----------------------------------------
145      !        EMP, SFX and QNS effects
146      !----------------------------------------
147      !                                   !==  Set before sbc tracer content fields  ==!
148      !
149      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
150         !
151         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
152              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
153            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
154            zfact = 0.5_wp
155            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
156            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
157         ELSE                                         ! No restart or restart not found: Euler forward time stepping
158            zfact = 1._wp
159            sbc_tsc_b(:,:,:) = 0._wp
160         ENDIF
161      ELSE                                         ! Swap of forcing fields
162         zfact = 0.5_wp
163         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
164      ENDIF
165      !                                          Compute now sbc tracer content fields
166      !                                          *************************************
167
168                                                   ! Concentration dilution effect on (t,s) due to 
169                                                   ! evaporation, precipitation and qns, but not river runoff
170                                               
171      IF( .NOT.ln_linssh ) THEN                    !* Variable Volume case  ==>> heat content of mass flux is in qns
172         DO jj = 1, jpj
173            DO ji = 1, jpi 
174               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
175               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
176            END DO
177         END DO
178      ELSE                                         !* Constant Volume case ==>> Concentration dilution effect
179         DO jj = 2, jpj
180            DO ji = fs_2, fs_jpim1   ! vector opt.
181               ! temperature : heat flux
182               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
183                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
184               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
185               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
186                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
187            END DO
188         END DO
189         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
190         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
191      ENDIF
192      !
193      DO jn = 1, jpts
194         DO jj = 2, jpj
195            DO ji = fs_2, fs_jpim1   ! vector opt.
196               z1_e3t = zfact / e3t_n(ji,jj,1)
197               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
198            END DO
199         END DO
200      END DO
201      !                                          Write in the ocean restart file
202      !                                          *******************************
203      IF( lrst_oce ) THEN
204         IF(lwp) WRITE(numout,*)
205         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
206            &                    'at it= ', kt,' date= ', ndastp
207         IF(lwp) WRITE(numout,*) '~~~~'
208         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
209         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
210      ENDIF
211      !
212      !
213      !----------------------------------------
214      !       Ice Shelf effects (ISF)
215      !     tbl treated as in Losh (2008) JGR
216      !----------------------------------------
217      !
218      IF( nn_isf > 0 ) THEN
219         zfact = 0.5_wp
220         DO jj = 2, jpj
221            DO ji = fs_2, fs_jpim1
222         
223               ikt = misfkt(ji,jj)
224               ikb = misfkb(ji,jj)
225   
226               ! level fully include in the ice shelf boundary layer
227               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
228               ! sign - because fwf sign of evapo (rnf sign of precip)
229               DO jk = ikt, ikb - 1
230               ! compute tfreez for the temperature correction (we add water at freezing temperature)
231               ! compute trend
232                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
233                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj)
234                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
235                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
236               END DO
237   
238               ! level partially include in ice shelf boundary layer
239               ! compute tfreez for the temperature correction (we add water at freezing temperature)
240               ! compute trend
241               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
242                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
243               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
244                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
245            END DO
246         END DO
247         IF( lrst_oce ) THEN
248            IF(lwp) WRITE(numout,*)
249            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
250               &                    'at it= ', kt,' date= ', ndastp
251            IF(lwp) WRITE(numout,*) '~~~~'
252            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
253            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
254            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
255         ENDIF
256      END IF
257      !
258      !----------------------------------------
259      !        River Runoff effects
260      !----------------------------------------
261      !
262      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
263         zfact = 0.5_wp
264         DO jj = 2, jpj 
265            DO ji = fs_2, fs_jpim1
266               IF( rnf(ji,jj) /= 0._wp ) THEN
267                  zdep = zfact / h_rnf(ji,jj)
268                  DO jk = 1, nk_rnf(ji,jj)
269                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
270                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
271                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
272                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
273                  END DO
274               ENDIF
275            END DO 
276         END DO 
277      ENDIF
278 
279      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
280         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
281         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
282         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
283         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
284         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
285      ENDIF
286      !
287      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
288         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
289      !
290      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
291      !
292   END SUBROUTINE tra_sbc
293
294   !!======================================================================
295END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.