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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 3952

Last change on this file since 3952 was 3952, checked in by flavoni, 11 years ago

ok for LIM2, still reproducibility error on LIM3 case

  • Property svn:keywords set to Id
File size: 13.0 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   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   tra_sbc      : update the tracer trend at ocean surface
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE dom_oce         ! ocean space domain variables
19   USE phycst          ! physical constant
20   USE traqsr          ! solar radiation penetration
21   USE trdmod_oce      ! ocean trends
22   USE trdtra          ! ocean trends
23   USE in_out_manager  ! I/O manager
24   USE prtctl          ! Print control
25   USE restart         ! ocean restart
26   USE sbcrnf          ! River runoff 
27   USE sbcmod          ! ln_rnf 
28   USE iom
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE wrk_nemo        ! Memory Allocation
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_sbc    ! routine called by step.F90
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE tra_sbc ( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE tra_sbc  ***
51      !!                   
52      !! ** Purpose :   Compute the tracer surface boundary condition trend of
53      !!      (flux through the interface, concentration/dilution effect)
54      !!      and add it to the general trend of tracer equations.
55      !!
56      !! ** Method :
57      !!      Following Roullet and Madec (2000), the air-sea flux can be divided
58      !!      into three effects: (1) Fext, external forcing;
59      !!      (2) Fwi, concentration/dilution effect due to water exchanged
60      !!         at the surface by evaporation, precipitations and runoff (E-P-R);
61      !!      (3) Fwe, tracer carried with the water that is exchanged.
62      !!
63      !!      Fext, flux through the air-sea interface for temperature and salt:
64      !!            - temperature : heat flux q (w/m2). If penetrative solar
65      !!         radiation q is only the non solar part of the heat flux, the
66      !!         solar part is added in traqsr.F routine.
67      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
68      !!            - salinity    : no salt flux
69      !!
70      !!      The formulation for Fwb and Fwi vary according to the free
71      !!      surface formulation (linear or variable volume).
72      !!      * Linear free surface
73      !!            The surface freshwater flux modifies the ocean volume
74      !!         and thus the concentration of a tracer and the temperature.
75      !!         First order of the effect of surface freshwater exchange
76      !!         for salinity, it can be neglected on temperature (especially
77      !!         as the temperature of precipitations and runoffs is usually
78      !!         unknown).
79      !!            - temperature : we assume that the temperature of both
80      !!         precipitations and runoffs is equal to the SST, thus there
81      !!         is no additional flux since in this case, the concentration
82      !!         dilution effect is balanced by the net heat flux associated
83      !!         to the freshwater exchange (Fwe+Fwi=0):
84      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
85      !!            - salinity    : evaporation, precipitation and runoff
86      !!         water has a zero salinity (Fwe=0), thus only Fwi remains:
87      !!            sa = sa + emp * sn / e3t   for k=1
88      !!         where emp, the surface freshwater budget (evaporation minus
89      !!         precipitation minus runoff) given in kg/m2/s is divided
90      !!         by 1035 kg/m3 (density of ocena water) to obtain m/s.   
91      !!         Note: even though Fwe does not appear explicitly for
92      !!         temperature in this routine, the heat carried by the water
93      !!         exchanged through the surface is part of the total heat flux
94      !!         forcing and must be taken into account in the global heat
95      !!         balance).
96      !!      * nonlinear free surface (variable volume, lk_vvl)
97      !!         contrary to the linear free surface case, Fwi is properly
98      !!         taken into account by using the true layer thicknesses to       
99      !!         calculate tracer content and advection. There is no need to
100      !!         deal with it in this routine.
101      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
102      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
103      !!
104      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
105      !!                with the tracer surface boundary condition
106      !!              - save the trend it in ttrd ('key_trdtra')
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
109      !!
110      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
111      REAL(wp) ::   zfact, z1_e3t, zsrau, zdep
112      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
113      !!----------------------------------------------------------------------
114      !
115      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
120         IF(lwp) WRITE(numout,*) '~~~~~~~ '
121      ENDIF
122
123      zsrau = 1. / rau0             ! initialization
124
125      IF( l_trdtra )   THEN                    !* Save ta and sa trends
126         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
127         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
128         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
129      ENDIF
130
131!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
132      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
133         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
134         qsr(:,:) = 0.e0                     ! qsr set to zero
135      ENDIF
136
137      !----------------------------------------
138      !        EMP, EMPS and QNS effects
139      !----------------------------------------
140      !                                          Set before sbc tracer content fields
141      !                                          ************************************
142      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
143         !                                         ! -----------------------------------
144         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
145              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
146            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
147            zfact = 0.5e0
148            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
149            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
150         ELSE                                         ! No restart or restart not found: Euler forward time stepping
151            zfact = 1.e0
152            sbc_tsc_b(:,:,:) = 0.e0
153         ENDIF
154      ELSE                                         ! Swap of forcing fields
155         !                                         ! ----------------------
156         zfact = 0.5e0
157         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
158      ENDIF
159      !                                          Compute now sbc tracer content fields
160      !                                          *************************************
161
162                                                   ! Concentration dilution effect on (t,s) due to 
163                                                   ! evaporation, precipitation and qns, but not river runoff
164                                               
165      IF( lk_vvl ) THEN                            ! Variable Volume case
166         DO jj = 1, jpj
167            DO ji = 1, jpi 
168               ! temperature : heat flux + cooling/heating effet of EMP flux
169               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
170               ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration
171#if defined key_lim3
172               sbc_tsc(ji,jj,jp_sal) = emps(ji,jj) * zsrau  ! IOVINO + CLEM
173#else
174               sbc_tsc(ji,jj,jp_sal) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal)
175#endif
176            END DO
177         END DO
178      ELSE                                         ! Constant Volume case
179         DO jj = 2, jpj
180            DO ji = fs_2, fs_jpim1   ! vector opt.
181               ! temperature : heat flux
182               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj)
183               ! salinity    : salt flux + concent./dilut. effect (both in emps)
184#if defined key_lim3
185               sbc_tsc(ji,jj,jp_sal) = zsrau * emps(ji,jj) + zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_sal) ! IOVINO
186#else
187               sbc_tsc(ji,jj,jp_sal) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal) ! regular expression
188#endif
189            END DO
190         END DO
191      ENDIF
192      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff 
193      DO jn = 1, jpts
194         DO jj = 2, jpj
195            DO ji = fs_2, fs_jpim1   ! vector opt.
196               z1_e3t = zfact / fse3t(ji,jj,1)
197               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
198            END DO
199         END DO
200      END DO
201      !                                          Write in the ocean restart file
202      !                                          *******************************
203      IF( lrst_oce ) THEN
204         IF(lwp) WRITE(numout,*)
205         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
206            &                    'at it= ', kt,' date= ', ndastp
207         IF(lwp) WRITE(numout,*) '~~~~'
208         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
209         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
210      ENDIF
211      !
212      !----------------------------------------
213      !        River Runoff effects
214      !----------------------------------------
215      !
216      zfact = 0.5e0
217
218      ! Effect on (t,s) due to river runoff (dilution effect automatically applied via vertical tracer advection)
219      IF( ln_rnf ) THEN 
220         DO jj = 2, jpj 
221            DO ji = fs_2, fs_jpim1
222               zdep = 1. / h_rnf(ji,jj)
223               zdep = zfact * zdep 
224               IF ( rnf(ji,jj) /= 0._wp ) THEN
225                  DO jk = 1, nk_rnf(ji,jj)
226                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
227                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
228                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
229                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
230                  END DO
231               ENDIF
232            END DO 
233         END DO 
234      ENDIF 
235!!gm  It should be useless
236      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )    ;    CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
237
238      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
239         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
240         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
241         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_nsr, ztrdt )
242         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_nsr, ztrds )
243         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
244      ENDIF
245      !
246      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
247         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
248      !
249      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
250      !
251   END SUBROUTINE tra_sbc
252
253   !!======================================================================
254END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.