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/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 5389

Last change on this file since 5389 was 5389, checked in by deazer, 9 years ago

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