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

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 5 years ago

Full set of changes as in the original branch

File size: 17.0 KB
RevLine 
[3]1MODULE trasbc
2   !!==============================================================================
3   !!                       ***  MODULE  trasbc  ***
4   !! Ocean active tracers:  surface boundary condition
5   !!==============================================================================
[2528]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
[5120]11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing
[3]12   !!----------------------------------------------------------------------
[503]13
14   !!----------------------------------------------------------------------
[3]15   !!   tra_sbc      : update the tracer trend at ocean surface
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and active tracers
[888]18   USE sbc_oce         ! surface boundary condition: ocean
[3]19   USE dom_oce         ! ocean space domain variables
20   USE phycst          ! physical constant
[4990]21   USE sbcmod          ! ln_rnf 
22   USE sbcrnf          ! River runoff 
[5431]23   USE sbcisf          ! Ice shelf   
[216]24   USE traqsr          ! solar radiation penetration
[4990]25   USE trd_oce         ! trends: ocean variables
26   USE trdtra          ! trends manager: tracers
[11134]27   USE tradwl          ! solar radiation penetration (downwell method)
[4990]28   !
[3]29   USE in_out_manager  ! I/O manager
[258]30   USE prtctl          ! Print control
[2528]31   USE iom
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[3294]33   USE wrk_nemo        ! Memory Allocation
34   USE timing          ! Timing
[4990]35   USE eosbn2
[11134]36#if defined key_asminc   
37   USE asminc          ! Assimilation increment
38#endif
[3]39
40   IMPLICIT NONE
41   PRIVATE
42
[503]43   PUBLIC   tra_sbc    ! routine called by step.F90
[3]44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
[4990]49   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
[888]50   !! $Id$
[2715]51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE tra_sbc ( kt )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE tra_sbc  ***
58      !!                   
59      !! ** Purpose :   Compute the tracer surface boundary condition trend of
60      !!      (flux through the interface, concentration/dilution effect)
61      !!      and add it to the general trend of tracer equations.
62      !!
63      !! ** Method :
[664]64      !!      Following Roullet and Madec (2000), the air-sea flux can be divided
65      !!      into three effects: (1) Fext, external forcing;
66      !!      (2) Fwi, concentration/dilution effect due to water exchanged
67      !!         at the surface by evaporation, precipitations and runoff (E-P-R);
68      !!      (3) Fwe, tracer carried with the water that is exchanged.
[3625]69      !!            - salinity    : salt flux only due to freezing/melting
70      !!            sa = sa +  sfx / rau0 / e3t  for k=1
[664]71      !!
72      !!      Fext, flux through the air-sea interface for temperature and salt:
[3]73      !!            - temperature : heat flux q (w/m2). If penetrative solar
74      !!         radiation q is only the non solar part of the heat flux, the
75      !!         solar part is added in traqsr.F routine.
76      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
77      !!            - salinity    : no salt flux
[664]78      !!
79      !!      The formulation for Fwb and Fwi vary according to the free
80      !!      surface formulation (linear or variable volume).
81      !!      * Linear free surface
82      !!            The surface freshwater flux modifies the ocean volume
[3]83      !!         and thus the concentration of a tracer and the temperature.
84      !!         First order of the effect of surface freshwater exchange
85      !!         for salinity, it can be neglected on temperature (especially
[664]86      !!         as the temperature of precipitations and runoffs is usually
87      !!         unknown).
[3]88      !!            - temperature : we assume that the temperature of both
89      !!         precipitations and runoffs is equal to the SST, thus there
90      !!         is no additional flux since in this case, the concentration
91      !!         dilution effect is balanced by the net heat flux associated
[664]92      !!         to the freshwater exchange (Fwe+Fwi=0):
93      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
[3]94      !!            - salinity    : evaporation, precipitation and runoff
[3625]95      !!         water has a zero salinity  but there is a salt flux due to
96      !!         freezing/melting, thus:
97      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
98      !!                    + sfx    / rau0 / e3t
[3]99      !!         where emp, the surface freshwater budget (evaporation minus
100      !!         precipitation minus runoff) given in kg/m2/s is divided
[4990]101      !!         by rau0 (density of sea water) to obtain m/s.   
[664]102      !!         Note: even though Fwe does not appear explicitly for
103      !!         temperature in this routine, the heat carried by the water
104      !!         exchanged through the surface is part of the total heat flux
105      !!         forcing and must be taken into account in the global heat
106      !!         balance).
107      !!      * nonlinear free surface (variable volume, lk_vvl)
108      !!         contrary to the linear free surface case, Fwi is properly
109      !!         taken into account by using the true layer thicknesses to       
110      !!         calculate tracer content and advection. There is no need to
111      !!         deal with it in this routine.
112      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
113      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
[3]114      !!
115      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
116      !!                with the tracer surface boundary condition
[4990]117      !!              - send trends to trdtra module (l_trdtra=T)
[503]118      !!----------------------------------------------------------------------
[2528]119      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[3]120      !!
[2528]121      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
[4990]122      INTEGER  ::   ikt, ikb 
123      INTEGER  ::   nk_isf
[3625]124      REAL(wp) ::   zfact, z1_e3t, zdep
[4990]125      REAL(wp) ::   zalpha, zhk
[3294]126      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
[3]127      !!----------------------------------------------------------------------
[3294]128      !
129      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
130      !
[3]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
[4990]137      IF( l_trdtra ) THEN                    !* Save ta and sa trends
[3294]138         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
139         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
140         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]141      ENDIF
142
[2528]143!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
[11134]144      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration
[2528]145         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
146         qsr(:,:) = 0.e0                     ! qsr set to zero
147      ENDIF
[3]148
[2528]149      !----------------------------------------
[4990]150      !        EMP, SFX and QNS effects
[2528]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'
[4990]159            zfact = 0.5_wp
[2528]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
[4990]163            zfact = 1._wp
164            sbc_tsc_b(:,:,:) = 0._wp
[2528]165         ENDIF
166      ELSE                                         ! Swap of forcing fields
167         !                                         ! ----------------------
[4990]168         zfact = 0.5_wp
[2528]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                                               
[3625]177      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
[3294]178         DO jj = 1, jpj
179            DO ji = 1, jpi 
[3625]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
[2528]182            END DO
[3]183         END DO
[3625]184      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
[2528]185         DO jj = 2, jpj
186            DO ji = fs_2, fs_jpim1   ! vector opt.
187               ! temperature : heat flux
[3625]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
[2528]193            END DO
194         END DO
[4990]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
[2528]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
[3]206      END DO
[2528]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      !
[4990]218      !
[2528]219      !----------------------------------------
[4990]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               ! compute trend
238                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
[5628]239                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj)
[4990]240                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
241                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
242               END DO
243   
244               ! level partially include in ice shelf boundary layer
245               ! compute tfreez for the temperature correction (we add water at freezing temperature)
246               ! compute trend
247               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
[5628]248                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
[4990]249               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
250                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
251            END DO
252         END DO
253         IF( lrst_oce ) THEN
254            IF(lwp) WRITE(numout,*)
255            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
256               &                    'at it= ', kt,' date= ', ndastp
257            IF(lwp) WRITE(numout,*) '~~~~'
258            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
259            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
260            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
261         ENDIF
262      END IF
263      !
264      !----------------------------------------
[2528]265      !        River Runoff effects
266      !----------------------------------------
267      !
[3764]268      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
269         zfact = 0.5_wp
[2528]270         DO jj = 2, jpj 
271            DO ji = fs_2, fs_jpim1
[3764]272               IF( rnf(ji,jj) /= 0._wp ) THEN
273                  zdep = zfact / h_rnf(ji,jj)
[2528]274                  DO jk = 1, nk_rnf(ji,jj)
275                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
276                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
277                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
278                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
[2715]279                  END DO
[2528]280               ENDIF
[2715]281            END DO 
282         END DO 
[3764]283      ENDIF
[11134]284
285#if defined key_asminc
286! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
287! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
288! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
289      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
290         DO jj = 2, jpj 
291            DO ji = fs_2, fs_jpim1
292               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
293               DO jk = 1, jpkm1
294                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
295                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
296                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
297                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
298               END DO
299            END DO 
300         END DO 
301      ENDIF
302#endif
[3764]303 
[4990]304      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
[2528]305         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
306         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[4990]307         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
308         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
[3294]309         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
[216]310      ENDIF
[503]311      !
[2528]312      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
313         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
[503]314      !
[3294]315      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
316      !
[3]317   END SUBROUTINE tra_sbc
318
319   !!======================================================================
320END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.