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

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 @ 8892

Last change on this file since 8892 was 8892, checked in by frrh, 6 years ago

Commit updates with debugging write statements.
File size: 13.8 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 eosbn2         ! Equation Of State
22   USE sbcmod         ! ln_rnf 
23   USE sbcrnf         ! River runoff 
24   USE sbcisf         ! Ice shelf   
25   USE iscplini       ! Ice sheet coupling
26   USE traqsr         ! solar radiation penetration
27   USE trd_oce        ! trends: ocean variables
28   USE trdtra         ! trends manager: tracers
29use lib_fortran
30   !
31   USE in_out_manager ! I/O manager
32   USE prtctl         ! Print control
33   USE iom            ! xIOS server
34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
35   USE wrk_nemo       ! Memory Allocation
36   USE timing         ! Timing
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   tra_sbc   ! routine called by step.F90
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_sbc ( kt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_sbc  ***
55      !!                   
56      !! ** Purpose :   Compute the tracer surface boundary condition trend of
57      !!      (flux through the interface, concentration/dilution effect)
58      !!      and add it to the general trend of tracer equations.
59      !!
60      !! ** Method :   The (air+ice)-sea flux has two components:
61      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface);
62      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.
63      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe,
64      !!             they are simply added to the tracer trend (tsa).
65      !!               In linear free surface case (ln_linssh=T), the volume of the
66      !!             ocean does not change with the water exchanges at the (air+ice)-sea
67      !!             interface. Therefore another term has to be added, to mimic the
68      !!             concentration/dilution effect associated with water exchanges.
69      !!
70      !! ** Action  : - Update tsa with the surface boundary condition trend
71      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
72      !!----------------------------------------------------------------------
73      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
74      !
75      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
76      INTEGER  ::   ikt, ikb              ! local integers
77      REAL(wp) ::   zfact, z1_e3t, zdep   ! local scalar
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
79      !!----------------------------------------------------------------------
80      !
81      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
82      !
83      IF( kt == nit000 ) THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
86         IF(lwp) WRITE(numout,*) '~~~~~~~ '
87      ENDIF
88      !
89      IF( l_trdtra ) THEN                    !* Save ta and sa trends
90         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
91         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
92         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
93      ENDIF
94
95write(numout,*) "RSRH tra_sbc tsb, tsa, tsn AA", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
96
97      !
98!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
99      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
100         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
101         qsr(:,:) = 0._wp                     ! qsr set to zero
102      ENDIF
103
104      !----------------------------------------
105      !        EMP, SFX and QNS effects
106      !----------------------------------------
107      !                             !==  Set before sbc tracer content fields  ==!
108      IF( kt == nit000 ) THEN             !* 1st time-step
109         IF( ln_rstart .AND.    &               ! Restart: read in restart file
110              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
111            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
112            zfact = 0.5_wp
113            sbc_tsc(:,:,:) = 0._wp
114            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
115            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
116         ELSE                                   ! No restart or restart not found: Euler forward time stepping
117            zfact = 1._wp
118            sbc_tsc(:,:,:) = 0._wp
119            sbc_tsc_b(:,:,:) = 0._wp
120         ENDIF
121      ELSE                                !* other time-steps: swap of forcing fields
122         zfact = 0.5_wp
123         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
124      ENDIF
125write(numout,*) "RSRH tra_sbc tsb, tsa, tsn BB", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
126      !                             !==  Now sbc tracer content fields  ==!
127      DO jj = 2, jpj
128         DO ji = fs_2, fs_jpim1   ! vector opt.
129            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
130            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
131         END DO
132      END DO
133      IF( ln_linssh ) THEN                !* linear free surface 
134         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem)
137               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal)
138            END DO
139         END DO                                 !==>> output c./d. term
140         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )
141         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )
142      ENDIF
143write(numout,*) "RSRH tra_sbc tsb, tsa, tsn CC", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
144write(numout,*) "RSRH tra_sbc qns, emp, sfx FF", kt, glob_sum(qns(:,:)), glob_sum(emp(:,:)), glob_sum(sfx(:,:)); flush(numout)
145write(numout,*) "RSRH tra_sbc tsn FF", kt, glob_sum(tsn(:,:,1,1)); flush(numout)
146write(numout,*) "RSRH tra_sbc e3t_n FF", kt, glob_sum(e3t_n(:,:,1)); flush(numout)
147write(numout,*) "RSRH tra_sbc sbc_tsc_b, sbc_tsc FF", kt, glob_sum(sbc_tsc_b(:,:,1)), glob_sum(sbc_tsc(:,:,1)); flush(numout)
148      !
149      DO jn = 1, jpts               !==  update tracer trend  ==!
150         DO jj = 2, jpj
151            DO ji = fs_2, fs_jpim1   ! vector opt. 
152               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1)
153            END DO
154         END DO
155      END DO
156write(numout,*) "RSRH tra_sbc tsb, tsa, tsn DD", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
157      !                 
158      IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==!
159         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
160         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
161      ENDIF
162      !
163      !----------------------------------------
164      !       Ice Shelf effects (ISF)
165      !     tbl treated as in Losh (2008) JGR
166      !----------------------------------------
167      !
168!!gm BUG ?   Why no differences between non-linear and linear free surface ?
169!!gm         probably taken into account in r1_hisf_tbl : to be verified
170      IF( ln_isf ) THEN
171write(numout,*) "RSRH tra_sbc risf_tsc_b, risf_tsc FF", kt, glob_sum(risf_tsc_b(:,:,1)), glob_sum(risf_tsc(:,:,1)); flush(numout)
172         zfact = 0.5_wp
173         DO jj = 2, jpj
174            DO ji = fs_2, fs_jpim1
175               !
176               ikt = misfkt(ji,jj)
177               ikb = misfkb(ji,jj)
178               !
179               ! level fully include in the ice shelf boundary layer
180               ! sign - because fwf sign of evapo (rnf sign of precip)
181               DO jk = ikt, ikb - 1
182               ! compute trend
183                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                &
184                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
185                     &           * r1_hisf_tbl(ji,jj)
186               END DO
187   
188               ! level partially include in ice shelf boundary layer
189               ! compute trend
190               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 &
191                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
192                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
193
194            END DO
195         END DO
196      END IF
197write(numout,*) "RSRH tra_sbc tsb, tsa, tsn EE", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
198      !
199      !----------------------------------------
200      !        River Runoff effects
201      !----------------------------------------
202      !
203      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff
204write(numout,*) "RSRH tra_sbc rnf_tsc_b, rnf_tsc FF", kt, glob_sum(rnf_tsc_b(:,:,1)), glob_sum(rnf_tsc(:,:,1)); flush(numout)
205         zfact = 0.5_wp
206         DO jj = 2, jpj 
207            DO ji = fs_2, fs_jpim1
208               IF( rnf(ji,jj) /= 0._wp ) THEN
209                  zdep = zfact / h_rnf(ji,jj)
210                  DO jk = 1, nk_rnf(ji,jj)
211                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 &
212                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
213                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 &
214                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
215                  END DO
216               ENDIF
217            END DO 
218         END DO 
219      ENDIF
220write(numout,*) "RSRH tra_sbc tsb, tsa, tsn FF", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
221
222      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
223      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
224
225      !
226      !----------------------------------------
227      !        Ice Sheet coupling imbalance correction to have conservation
228      !----------------------------------------
229      !
230      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff
231         DO jk = 1,jpk
232            DO jj = 2, jpj 
233               DO ji = fs_2, fs_jpim1
234                  zdep = 1._wp / e3t_n(ji,jj,jk) 
235                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem)                       &
236                      &                                         * zdep
237                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal)                       &
238                      &                                         * zdep 
239               END DO 
240            END DO 
241         END DO
242      ENDIF
243write(numout,*) "RSRH tra_sbc tsb, tsa, tsn GG", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
244
245      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
246         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
247         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
248         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
249         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
250         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
251      ENDIF
252      !
253      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
254         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
255      !
256      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
257write(numout,*) "RSRH tra_sbc tsb, tsa, tsn HH", kt, glob_sum(tsb(:,:,:,1)), glob_sum(tsa(:,:,:,1)),glob_sum(tsn(:,:,:,1)); flush(numout)
258      !
259   END SUBROUTINE tra_sbc
260
261   !!======================================================================
262END MODULE trasbc
Note: See TracBrowser for help on using the repository browser.