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

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

File size: 17.7 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
27   !
[3]28   USE in_out_manager  ! I/O manager
[258]29   USE prtctl          ! Print control
[2528]30   USE iom
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[3294]32   USE wrk_nemo        ! Memory Allocation
33   USE timing          ! Timing
[4990]34   USE eosbn2
[8400]35#if defined key_asminc   
36   USE asminc          ! Assimilation increment
37#endif
[3]38
39   IMPLICIT NONE
40   PRIVATE
41
[503]42   PUBLIC   tra_sbc    ! routine called by step.F90
[3]43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
[4990]48   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
[888]49   !! $Id$
[2715]50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]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 :
[664]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.
[3625]68      !!            - salinity    : salt flux only due to freezing/melting
69      !!            sa = sa +  sfx / rau0 / e3t  for k=1
[664]70      !!
71      !!      Fext, flux through the air-sea interface for temperature and salt:
[3]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
[664]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
[3]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
[664]85      !!         as the temperature of precipitations and runoffs is usually
86      !!         unknown).
[3]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
[664]91      !!         to the freshwater exchange (Fwe+Fwi=0):
92      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
[3]93      !!            - salinity    : evaporation, precipitation and runoff
[3625]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
[3]98      !!         where emp, the surface freshwater budget (evaporation minus
99      !!         precipitation minus runoff) given in kg/m2/s is divided
[4990]100      !!         by rau0 (density of sea water) to obtain m/s.   
[664]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.
[3]113      !!
114      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
115      !!                with the tracer surface boundary condition
[4990]116      !!              - send trends to trdtra module (l_trdtra=T)
[503]117      !!----------------------------------------------------------------------
[2528]118      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[3]119      !!
[2528]120      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
[4990]121      INTEGER  ::   ikt, ikb 
122      INTEGER  ::   nk_isf
[3625]123      REAL(wp) ::   zfact, z1_e3t, zdep
[4990]124      REAL(wp) ::   zalpha, zhk
[8400]125      REAL(wp) ::  zt_frz, zpress
[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,*) '~~~~~~~ '
[10774]135         IF(lwp .AND. lflush) CALL flush(numout)
[3]136      ENDIF
137
[4990]138      IF( l_trdtra ) THEN                    !* Save ta and sa trends
[3294]139         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
140         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
141         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]142      ENDIF
143
[2528]144!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
145      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
146         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
147         qsr(:,:) = 0.e0                     ! qsr set to zero
148      ENDIF
[3]149
[2528]150      !----------------------------------------
[4990]151      !        EMP, SFX and QNS effects
[2528]152      !----------------------------------------
153      !                                          Set before sbc tracer content fields
154      !                                          ************************************
155      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
156         !                                         ! -----------------------------------
157         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
158              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
159            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
[10774]160            IF(lwp .AND. lflush) CALL flush(numout)
[4990]161            zfact = 0.5_wp
[2528]162            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
163            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
164         ELSE                                         ! No restart or restart not found: Euler forward time stepping
[4990]165            zfact = 1._wp
[6793]166            sbc_tsc(:,:,:) = 0._wp
[4990]167            sbc_tsc_b(:,:,:) = 0._wp
[2528]168         ENDIF
169      ELSE                                         ! Swap of forcing fields
170         !                                         ! ----------------------
[4990]171         zfact = 0.5_wp
[2528]172         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
173      ENDIF
174      !                                          Compute now sbc tracer content fields
175      !                                          *************************************
176
177                                                   ! Concentration dilution effect on (t,s) due to 
178                                                   ! evaporation, precipitation and qns, but not river runoff
179                                               
[3625]180      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
[3294]181         DO jj = 1, jpj
182            DO ji = 1, jpi 
[3625]183               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
184               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
[2528]185            END DO
[3]186         END DO
[3625]187      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
[2528]188         DO jj = 2, jpj
189            DO ji = fs_2, fs_jpim1   ! vector opt.
190               ! temperature : heat flux
[3625]191               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
192                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
193               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
194               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
195                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
[2528]196            END DO
197         END DO
[4990]198         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
199         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
[2528]200      ENDIF
201      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff 
202      DO jn = 1, jpts
203         DO jj = 2, jpj
204            DO ji = fs_2, fs_jpim1   ! vector opt.
205               z1_e3t = zfact / fse3t(ji,jj,1)
206               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
207            END DO
208         END DO
[3]209      END DO
[2528]210      !                                          Write in the ocean restart file
211      !                                          *******************************
212      IF( lrst_oce ) THEN
[10770]213         IF(lwp .AND. nprint > 0) THEN
214            WRITE(numout,*)
215            WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
[2528]216            &                    'at it= ', kt,' date= ', ndastp
[10770]217            WRITE(numout,*) '~~~~'
[10774]218            IF(lflush) CALL flush(numout)
[10770]219         ENDIF
[9321]220         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
[2528]221         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
222         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
[9321]223         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
[2528]224      ENDIF
225      !
[4990]226      !
[2528]227      !----------------------------------------
[4990]228      !       Ice Shelf effects (ISF)
229      !     tbl treated as in Losh (2008) JGR
230      !----------------------------------------
231      !
232      IF( nn_isf > 0 ) THEN
233         zfact = 0.5e0
234         DO jj = 2, jpj
235            DO ji = fs_2, fs_jpim1
236         
237               ikt = misfkt(ji,jj)
238               ikb = misfkb(ji,jj)
239   
240               ! level fully include in the ice shelf boundary layer
241               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
242               ! sign - because fwf sign of evapo (rnf sign of precip)
243               DO jk = ikt, ikb - 1
244               ! compute tfreez for the temperature correction (we add water at freezing temperature)
245               ! compute trend
246                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
[6487]247                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj)
[4990]248                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
249                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
250               END DO
251   
252               ! level partially include in ice shelf boundary layer
253               ! compute tfreez for the temperature correction (we add water at freezing temperature)
254               ! compute trend
255               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
[6487]256                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
[4990]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
[10770]262            IF(lwp .AND. nprint > 0) THEN
263                WRITE(numout,*)
264                WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
[4990]265               &                    'at it= ', kt,' date= ', ndastp
[10770]266                WRITE(numout,*) '~~~~'
[10774]267                IF(lflush) CALL flush(numout)
[10770]268            ENDIF
[9321]269            IF(nn_timing == 2)  CALL timing_start('iom_rstput')
[4990]270            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
271            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
272            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
[9321]273            IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
[4990]274         ENDIF
275      END IF
276      !
277      !----------------------------------------
[2528]278      !        River Runoff effects
279      !----------------------------------------
280      !
[3764]281      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
282         zfact = 0.5_wp
[2528]283         DO jj = 2, jpj 
284            DO ji = fs_2, fs_jpim1
[3764]285               IF( rnf(ji,jj) /= 0._wp ) THEN
286                  zdep = zfact / h_rnf(ji,jj)
[2528]287                  DO jk = 1, nk_rnf(ji,jj)
288                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
289                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
290                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
291                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
[2715]292                  END DO
[2528]293               ENDIF
[2715]294            END DO 
295         END DO 
[3764]296      ENDIF
[7993]297
298      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
299      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
300
[8400]301#if defined key_asminc
302! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
303! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
304! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
305      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
306         DO jj = 2, jpj 
307            DO ji = fs_2, fs_jpim1
308               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
309               DO jk = 1, jpkm1
310                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
311                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
312                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
313                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
314               END DO
315            END DO 
316         END DO 
317      ENDIF
318#endif
319 
[4990]320      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
[2528]321         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
322         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[4990]323         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
324         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
[3294]325         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
[216]326      ENDIF
[503]327      !
[2528]328      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
329         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
[503]330      !
[3294]331      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
332      !
[3]333   END SUBROUTINE tra_sbc
334
335   !!======================================================================
336END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.