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_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 6012

Last change on this file since 6012 was 6012, checked in by mathiot, 8 years ago

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