source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 5605

Last change on this file since 5605 was 5605, checked in by mathiot, 6 years ago

Marine glacier termini: initial commit

  • Property svn:keywords set to Id
File size: 16.7 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 traicw          ! ice wall parametrization
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      INTEGER  ::   nk_isf
121      REAL(wp) ::   zfact, z1_e3t, zdep
122      REAL(wp) ::   zalpha, zhk
123      REAL(wp) ::  zt_frz, zpress
124      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
125      !!----------------------------------------------------------------------
126      !
127      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
128      !
129      IF( kt == nit000 ) THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
132         IF(lwp) WRITE(numout,*) '~~~~~~~ '
133      ENDIF
134
135      IF( l_trdtra ) THEN                    !* Save ta and sa trends
136         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
137         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
138         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
139      ENDIF
140
141!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
142      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
143         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
144         qsr(:,:) = 0.e0                     ! qsr set to zero
145      ENDIF
146
147      !----------------------------------------
148      !        EMP, SFX and QNS effects
149      !----------------------------------------
150      !                                          Set before sbc tracer content fields
151      !                                          ************************************
152      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
153         !                                         ! -----------------------------------
154         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
155              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
156            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
157            zfact = 0.5_wp
158            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
159            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
160         ELSE                                         ! No restart or restart not found: Euler forward time stepping
161            zfact = 1._wp
162            sbc_tsc_b(:,:,:) = 0._wp
163         ENDIF
164      ELSE                                         ! Swap of forcing fields
165         !                                         ! ----------------------
166         zfact = 0.5_wp
167         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
168      ENDIF
169      !                                          Compute now sbc tracer content fields
170      !                                          *************************************
171
172                                                   ! Concentration dilution effect on (t,s) due to 
173                                                   ! evaporation, precipitation and qns, but not river runoff
174                                               
175      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
176         DO jj = 1, jpj
177            DO ji = 1, jpi 
178               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
179               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
180            END DO
181         END DO
182      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
183         DO jj = 2, jpj
184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               ! temperature : heat flux
186               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
187                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
188               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
189               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
190                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
191            END DO
192         END DO
193         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
194         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
195      ENDIF
196      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff 
197      DO jn = 1, jpts
198         DO jj = 2, jpj
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               z1_e3t = zfact / fse3t(ji,jj,1)
201               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
202            END DO
203         END DO
204      END DO
205      !                                          Write in the ocean restart file
206      !                                          *******************************
207      IF( lrst_oce ) THEN
208         IF(lwp) WRITE(numout,*)
209         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
210            &                    'at it= ', kt,' date= ', ndastp
211         IF(lwp) WRITE(numout,*) '~~~~'
212         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
213         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
214      ENDIF
215      !
216      !
217      !----------------------------------------
218      !       Ice Shelf effects (ISF)
219      !     tbl treated as in Losh (2008) JGR
220      !----------------------------------------
221      !
222      IF( nn_isf > 0 ) THEN
223         zfact = 0.5e0
224         DO jj = 2, jpj
225            DO ji = fs_2, fs_jpim1
226         
227               ikt = misfkt(ji,jj)
228               ikb = misfkb(ji,jj)
229   
230               ! level fully include in the ice shelf boundary layer
231               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
232               ! sign - because fwf sign of evapo (rnf sign of precip)
233               DO jk = ikt, ikb - 1
234               ! compute tfreez for the temperature correction (we add water at freezing temperature)
235!                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
236                  zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
237               ! compute trend
238                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
239                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
240                     &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &
241                     &           * r1_hisf_tbl(ji,jj)
242                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
243                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
244               END DO
245   
246               ! level partially include in ice shelf boundary layer
247               ! compute tfreez for the temperature correction (we add water at freezing temperature)
248!               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04
249               zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
250               ! compute trend
251               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
252                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
253                  &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
254                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
255               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
256                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
257            END DO
258         END DO
259         IF( lrst_oce ) THEN
260            IF(lwp) WRITE(numout,*)
261            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
262               &                    'at it= ', kt,' date= ', ndastp
263            IF(lwp) WRITE(numout,*) '~~~~'
264            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
265            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
266            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
267         ENDIF
268      END IF
269      !
270      !----------------------------------------
271      !        River Runoff effects
272      !----------------------------------------
273      !
274      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
275         zfact = 0.5_wp
276         DO jj = 2, jpj 
277            DO ji = fs_2, fs_jpim1
278               IF( rnf(ji,jj) /= 0._wp ) THEN
279                  zdep = zfact / h_rnf(ji,jj)
280                  DO jk = 1, nk_rnf(ji,jj)
281                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
282                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
283                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
284                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
285                  END DO
286               ENDIF
287            END DO 
288         END DO 
289      ENDIF
290
291
292      ! compute trend
293      IF ( ln_traicw ) THEN
294         tsa(:,:,:,1) = tsa(:,:,:,1) + 0.5_wp * ( ricw_tsc(:,:,:,1) + ricw_tsc_b(:,:,:,1) )
295         tsa(:,:,:,2) = tsa(:,:,:,2) + 0.5_wp * ( ricw_tsc(:,:,:,2) + ricw_tsc_b(:,:,:,2) )
296      END IF
297 
298      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
299         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
300         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
301         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
302         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
303         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
304      ENDIF
305      !
306      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
307         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
308      !
309      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
310      !
311   END SUBROUTINE tra_sbc
312
313   !!======================================================================
314END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.