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

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90 @ 2068

Last change on this file since 2068 was 2068, checked in by mlelod, 14 years ago

ticket: #663 ensuring restartability and conservation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.6 KB
Line 
1MODULE trasbc
2   !!==============================================================================
3   !!                       ***  MODULE  trasbc  ***
4   !! Ocean active tracers:  surface boundary condition
5   !!==============================================================================
6   !! History :  OPA  !  2998-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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   tra_sbc      : update the tracer trend at ocean surface
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and active tracers
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE dom_oce         ! ocean space domain variables
18   USE phycst          ! physical constant
19   USE traqsr          ! solar radiation penetration
20   USE trdmod          ! ocean trends
21   USE trdmod_oce      ! ocean variables trends
22   USE iom
23   USE in_out_manager  ! I/O manager
24   USE restart         ! ocean restart
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   tra_sbc    ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE tra_sbc ( kt )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE tra_sbc  ***
46      !!                   
47      !! ** Purpose :   Compute the tracer surface boundary condition trend of
48      !!      (flux through the interface, concentration/dilution effect)
49      !!      and add it to the general trend of tracer equations.
50      !!
51      !! ** Method  :   Following Roullet and Madec (2000), the air-sea flux
52      !!      can be divided into three effects:
53      !!      (1) Fext, external forcing;
54      !!      (2) Fwi, concentration/dilution effect due to water exchanged
55      !!         at the surface by evaporation, precipitations and runoff (E-P-R);
56      !!      (3) Fwe, tracer carried with the water that is exchanged.
57      !!
58      !!      Fext, flux through the air-sea interface for temperature and salt:
59      !!            - temperature : heat flux q (w/m2). If penetrative solar
60      !!         radiation q is only the non solar part of the heat flux, the
61      !!         solar part is added in traqsr.F routine.
62      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
63      !!            - salinity    : only salt flux exchanged with sea-ice
64      !!
65      !!      The formulation for Fwb and Fwi vary according to the free
66      !!      surface formulation (linear or variable volume).
67      !!      * Linear free surface
68      !!            The surface freshwater flux modifies the ocean volume
69      !!         and thus the concentration of a tracer and the temperature.
70      !!         First order of the effect of surface freshwater exchange
71      !!         for salinity, it can be neglected on temperature (especially
72      !!         as the temperature of precipitations and runoffs is usually
73      !!         unknown).
74      !!            - temperature : we assume that the temperature of both
75      !!         precipitations and runoffs is equal to the SST, thus there
76      !!         is no additional flux since in this case, the concentration
77      !!         dilution effect is balanced by the net heat flux associated
78      !!         to the freshwater exchange (Fwe+Fwi=0):
79      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
80      !!            - salinity    : evaporation, precipitation and runoff
81      !!         water has a zero salinity (Fwe=0), thus only Fwi remains:
82      !!            sa = sa + emp * sn / e3t   for k=1
83      !!         where emp, the surface freshwater budget (evaporation minus
84      !!         precipitation minus runoff) given in kg/m2/s is divided
85      !!         by 1035 kg/m3 (density of ocena water) to obtain m/s.   
86      !!         Note: even though Fwe does not appear explicitly for
87      !!         temperature in this routine, the heat carried by the water
88      !!         exchanged through the surface is part of the total heat flux
89      !!         forcing and must be taken into account in the global heat
90      !!         balance).
91      !!      * nonlinear free surface (variable volume, lk_vvl)
92      !!         contrary to the linear free surface case, Fwi is properly
93      !!         taken into account by using the true layer thicknesses to       
94      !!         calculate tracer content and advection. There is no need to
95      !!         deal with it in this routine.
96      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
97      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
98      !!
99      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
100      !!                with the tracer surface boundary condition
101      !!              - save the trend it in ttrd ('key_trdtra')
102      !!----------------------------------------------------------------------
103      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
104      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
105      !!
106      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
107      !!
108      INTEGER  ::   ji, jj            ! dummy loop indices
109      REAL(wp) ::   zsrau, zse3t      ! temporary scalars
110      !!----------------------------------------------------------------------
111
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
115         IF(lwp) WRITE(numout,*) '~~~~~~~ '
116      ENDIF
117
118      zsrau = 1. / rau0             ! initialization
119#if defined key_zco
120      zse3t = 1. / e3t_0(1)
121#endif
122
123      IF( l_trdtra ) THEN           ! Save ta and sa trends
124         ztrdt(:,:,:) = ta(:,:,:) 
125         ztrds(:,:,:) = sa(:,:,:) 
126      ENDIF
127
128!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
129      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
130         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
131         qsr(:,:) = 0.e0                     ! qsr set to zero
132         IF( kt == nit000 ) THEN             ! idem on before field at nit000
133            qsr_b(:,:) = 0.e0                     
134         ENDIF
135      ENDIF
136      !                                            ! ---------------------------------------- !
137      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
138         !                                         ! ---------------------------------------- !
139         sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:)                         ! Swap the ocean forcing fields except at nit000
140         IF ( .NOT. lk_vvl ) sbc_trd_sc_b(:,:)   = sbc_trd_sc_n(:,:)
141      ENDIF
142      !                                            ! ---------------------------------------- !
143      IF( kt == nit000 ) THEN                      !   set the forcing field at nit000 - 1    !
144         !                                         ! ---------------------------------------- !
145         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
146            & iom_varid( numror, 'sbc_trd_hc_b', ldstop = .FALSE. ) > 0 ) THEN
147            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
148            CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend
149            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux
150            IF ( .NOT. lk_vvl ) THEN
151               CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend
152            ENDIF
153         ENDIF
154      ENDIF
155      !                                            ! ---------------------- !
156      IF( lk_vvl ) THEN                            !  Variable Volume case  !
157         !                                         ! ---------------------- !
158!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to sea-ice freezing/melting
159!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.s
160         IF ( neuler == 0 .AND. kt == nit000 ) THEN
161            DO jj = 2, jpj
162               DO ji = fs_2, fs_jpim1   ! vector opt.
163                  ! temperature : heat flux + cooling/heating effet of EMP flux
164                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1)
165#if ! defined key_zco
166                  zse3t = 1. / fse3t(ji,jj,1)
167#endif
168                  ta(ji,jj,1) = ta(ji,jj,1) + zse3t * sbc_trd_hc_n(ji,jj)
169                END DO
170            END DO
171         ELSE
172            DO jj = 2, jpj
173               DO ji = fs_2, fs_jpim1   ! vector opt.
174                  ! temperature : heat flux + cooling/heating effet of EMP flux
175                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1)
176#if ! defined key_zco
177                  zse3t = 1. / fse3t(ji,jj,1)
178#endif
179                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t
180               END DO
181            END DO
182         ENDIF
183         !                                         ! ---------------------- !
184      ELSE                                         !  Constant Volume case  !
185         !                                         ! ---------------------- !
186         IF ( neuler == 0 .AND. kt == nit000 ) THEN
187            DO jj = 2, jpj
188               DO ji = fs_2, fs_jpim1   ! vector opt.
189                  ! temperature : heat flux
190                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj)
191                  ! salinity    : salt flux + concent./dilut. effect (both in emps)
192                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1)
193#if ! defined key_zco
194                  zse3t = 1. / fse3t(ji,jj,1)
195#endif
196                  ta(ji,jj,1) = ta(ji,jj,1) + sbc_trd_hc_n(ji,jj) * zse3t
197                  sa(ji,jj,1) = sa(ji,jj,1) + sbc_trd_sc_n(ji,jj) * zse3t
198               END DO
199            END DO
200         ELSE
201            DO jj = 2, jpj
202               DO ji = fs_2, fs_jpim1   ! vector opt.
203                  ! temperature : heat flux
204                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj)
205                  ! salinity    : salt flux + concent./dilut. effect (both in emps)
206                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1)
207#if ! defined key_zco
208                  zse3t = 1. / fse3t(ji,jj,1)
209#endif
210                  ! temperature : heat flux
211                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t 
212                  sa(ji,jj,1) = sa(ji,jj,1) + 0.5 * ( sbc_trd_sc_b(ji,jj) + sbc_trd_sc_n(ji,jj) ) * zse3t
213               END DO
214            END DO
215         ENDIF
216         !
217      ENDIF
218
219      !                                            ! ---------------------------------------- !
220      IF( lrst_oce ) THEN                          !      Write in the ocean restart file     !
221         !                                         ! ---------------------------------------- !
222         IF(lwp) WRITE(numout,*)
223         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
224            &                    'at it= ', kt,' date= ', ndastp
225         IF(lwp) WRITE(numout,*) '~~~~'
226         CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_n )
227         IF ( .NOT. lk_vvl ) CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_n )
228         !
229      ENDIF
230
231      IF( l_trdtra ) THEN           ! save the sbc trends for diagnostic
232         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
233         ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
234         CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt)
235      ENDIF
236      !                             ! control print
237      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   &
238         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
239   END SUBROUTINE tra_sbc
240
241   !!======================================================================
242END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.