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.
tranxt.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/TRA/tranxt.F90 @ 15608

Last change on this file since 15608 was 15315, checked in by hadjt, 3 years ago

tracnticing - not allowing gridboxes to cool below freezing.

This was an option within CO6, which I have added again.

Namelist added: nam_tranxticing with nn_tranxticing

The version is committed in restrospect, so may not include all the changes to work.
File size: 22.3 KB
Line 
1MODULE tranxt
2   !!======================================================================
3   !!                       ***  MODULE  tranxt  ***
4   !! Ocean active tracers:  time stepping on temperature and salinity
5   !!======================================================================
6   !! History :  OPA  !  1991-11  (G. Madec)  Original code
7   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
8   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
9   !!             -   !  1996-04  (A. Weaver)  Euler forward step
10   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
11   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
13   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
14   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
15   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
16   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
18   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_nxt       : time stepping on tracers
23   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
24   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers variables
27   USE dom_oce         ! ocean space and time domain variables
28   USE sbc_oce         ! surface boundary condition: ocean
29   USE sbcrnf          ! river runoffs
30   USE sbcisf          ! ice shelf melting
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE trd_oce         ! trends: ocean variables
34   USE trdtra          ! trends manager: tracers
35   USE traqsr          ! penetrative solar radiation (needed for nksr)
36   USE phycst          ! physical constant
37   USE ldftra          ! lateral physics : tracers
38   USE ldfslp          ! lateral physics : slopes
39   USE bdy_oce  , ONLY : ln_bdy
40   USE bdytra          ! open boundary condition (bdy_tra routine)
41   !
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE timing          ! Timing
46#if defined key_agrif
47   USE agrif_oce_interp
48#endif
49   !JT for ctrl warn
50   USE lib_mpp
51   !JT for ctrl warn
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC   tra_nxt       ! routine called by step.F90
57   PUBLIC   tra_nxt_fix   ! to be used in trcnxt
58   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
59
60   !JT
61   INTEGER  ::   warn_1, warn_2   ! indicators for warning statement
62   INTEGER , PUBLIC ::   nn_tranxticing  ! region mean calculation
63   !JT
64
65   !! * Substitutions
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
69   !! $Id$
70   !! Software governed by the CeCILL license (see ./LICENSE)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE tra_nxt( kt )
75      !!----------------------------------------------------------------------
76      !!                   ***  ROUTINE tranxt  ***
77      !!
78      !! ** Purpose :   Apply the boundary condition on the after temperature 
79      !!             and salinity fields, achieved the time stepping by adding
80      !!             the Asselin filter on now fields and swapping the fields.
81      !!
82      !! ** Method  :   At this stage of the computation, ta and sa are the
83      !!             after temperature and salinity as the time stepping has
84      !!             been performed in trazdf_imp or trazdf_exp module.
85      !!
86      !!              - Apply lateral boundary conditions on (ta,sa)
87      !!             at the local domain   boundaries through lbc_lnk call,
88      !!             at the one-way open boundaries (ln_bdy=T),
89      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
90      !!
91      !!              - Update lateral boundary conditions on AGRIF children
92      !!             domains (lk_agrif=T)
93      !!
94      !! ** Action  : - tsb & tsn ready for the next time step
95      !!----------------------------------------------------------------------
96      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
97      !!
98      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
99      REAL(wp) ::   zfact            ! local scalars
100      ! JT
101      INTEGER  ::   ios   ! dummy loop indices
102      REAL(wp) ::   zfreeze   ! local scalars
103      ! JT
104      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
105      !!----------------------------------------------------------------------
106      !
107      IF( ln_timing )   CALL timing_start( 'tra_nxt')
108      !
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
112         IF(lwp) WRITE(numout,*) '~~~~~~~'
113
114
115         
116
117          NAMELIST/nam_tranxticing/ nn_tranxticing
118
119          REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
120          READ   ( numnam_ref, nam_tranxticing, IOSTAT=ios, ERR= 901 )
121901       IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tranxticing in reference namelist' )
122
123          REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
124          READ  ( numnam_cfg, nam_tranxticing, IOSTAT = ios, ERR = 902 )
125902       IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tranxticing in configuration namelist' )
126          IF(lwm) WRITE ( numond, nam_tranxticing )
127
128          IF(lwp) THEN                   ! Control print
129              WRITE(numout,*)
130              WRITE(numout,*) 'tra_nxt : Catching Icing'
131              WRITE(numout,*) '~~~~~~~~~~~~'
132              WRITE(numout,*) 'Namelist nam_tranxticing : catch freezing points '
133              WRITE(numout,*) 'Int switch for icing catch (0=Off; 1=Error; 2 Warn (verbose); 3 Warn (Quiet); nn_tranxticing  = ', nn_tranxticing
134          ENDIF
135     
136         
137 
138
139
140      ENDIF
141
142      ! Update after tracer on domain lateral boundaries
143      !
144#if defined key_agrif
145      CALL Agrif_tra                     ! AGRIF zoom boundaries
146#endif
147      !                                              ! local domain boundaries  (T-point, unchanged sign)
148      CALL lbc_lnk_multi( 'tranxt', tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1. )
149      !
150      IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
151 
152      ! set time step size (Euler/Leapfrog)
153      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
154      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
155      ENDIF
156
157      ! trends computation initialisation
158      IF( l_trdtra )   THEN                   
159         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
160         ztrdt(:,:,jpk) = 0._wp
161         ztrds(:,:,jpk) = 0._wp
162         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
163            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
164            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
165         ENDIF
166         ! total trend for the non-time-filtered variables.
167         zfact = 1.0 / rdt
168         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
169         DO jk = 1, jpkm1
170            ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
171            ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) / e3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
172         END DO
173         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
174         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
175         IF( ln_linssh ) THEN       ! linear sea surface height only
176            ! Store now fields before applying the Asselin filter
177            ! in order to calculate Asselin filter trend later.
178            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
179            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
180         ENDIF
181      ENDIF
182
183      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
184         DO jn = 1, jpts
185            DO jk = 1, jpkm1
186               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)   
187            END DO
188         END DO
189         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
190            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
191            ztrdt(:,:,:) = 0._wp
192            ztrds(:,:,:) = 0._wp
193            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
194            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
195         END IF
196         !
197      ELSE                                            ! Leap-Frog + Asselin filter time stepping
198         !
199         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface
200         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
201           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
202         ENDIF
203         !
204         CALL lbc_lnk_multi( 'tranxt', tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., &
205                  &          tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., &
206                  &          tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  )
207         !
208      ENDIF     
209      !JT
210      !
211#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice )
212
213      IF (nn_tranxticing .GT. 0) THEN
214        IF ( kt == nit000 ) warn_1=0 
215        warn_2=0 
216        DO jk = 1, jpkm1 
217          DO jj = 1, jpj 
218            DO ji = 1, jpi 
219              IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN 
220                ! calculate the freezing point
221                zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt (Abs(tsn(ji,jj,jk,jp_sal)))   & 
222                          - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + gdept_n(ji,jj,jk) ) 
223                IF ( tsn(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN
224                    IF (nn_tranxticing .lt. 3) THEN
225                 
226                      WRITE(numout,300) ' tranxticing: ', & 
227                      &   kt,zfreeze,tsn(ji,jj,jk,jp_tem),ji,jj,jk,narea,ji-1+njmpp-1+jpjglo,jj-1+nimpp-1+jpiglo
228                      !JT &   kt,zfreeze,tsn(ji,jj,jk,jp_tem),ji,jj,jk,narea,ji-1+njmpp-1+jpjzoom,jj-1+nimpp-1+jpizoom
229
230300 FORMAT(A14,1X,I7,1X,f9.2,1X,f9.2,1X,I4,1X,I4,1X,I4,1X,I4,1X,I4,1X,I4)
231                                       
232                    ENDIF
233                    tsn(ji,jj,jk,jp_tem)=zfreeze 
234                    warn_2=1 
235                ENDIF
236              ENDIF
237            END DO
238          END DO
239        END DO
240
241
242      !CALL mpp_max('warn_1')
243      !CALL mpp_max('warn_2')
244      !  IF ( (warn_1 == 0) .AND. (warn_2 /= 0) ) THEN
245      !    IF(lwp) THEN
246      !   
247      !      IF (nn_tranxticing .GT. 1) CALL ctl_warn( ' tranxticing: Temperatures dropping below freezing point, ', &
248      !                &    ' being forced to freezing point, no longer conservative' )
249      !               
250      !      IF (nn_tranxticing .EQ. 1) CALL ctl_stop( ' tranxticing: Temperatures dropping below freezing point, ', &
251      !                &    ' being forced to freezing point, no longer conservative' )
252      !               
253      !    ENDIF
254      !    warn_1=1
255      !  ENDIF
256
257      ENDIF 
258#endif 
259      !JT
260     
261      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
262         zfact = 1._wp / r2dt             
263         DO jk = 1, jpkm1
264            ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
265            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
266         END DO
267         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
268         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
269      END IF
270      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
271      !
272      !                        ! control print
273      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
274         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
275      !
276      IF( ln_timing )   CALL timing_stop('tra_nxt')
277      !
278   END SUBROUTINE tra_nxt
279
280
281   SUBROUTINE tra_nxt_fix( kt, kit000, cdtype, ptb, ptn, pta, kjpt )
282      !!----------------------------------------------------------------------
283      !!                   ***  ROUTINE tra_nxt_fix  ***
284      !!
285      !! ** Purpose :   fixed volume: apply the Asselin time filter and
286      !!                swap the tracer fields.
287      !!
288      !! ** Method  : - Apply a Asselin time filter on now fields.
289      !!              - swap tracer fields to prepare the next time_step.
290      !!
291      !! ** Action  : - tsb & tsn ready for the next time step
292      !!----------------------------------------------------------------------
293      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
294      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
295      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
296      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
297      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
298      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
299      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
300      !
301      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
302      REAL(wp) ::   ztn, ztd         ! local scalars
303      !!----------------------------------------------------------------------
304      !
305      IF( kt == kit000 )  THEN
306         IF(lwp) WRITE(numout,*)
307         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
308         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
309      ENDIF
310      !
311      DO jn = 1, kjpt
312         !
313         DO jk = 1, jpkm1
314            DO jj = 2, jpjm1
315               DO ji = fs_2, fs_jpim1
316                  ztn = ptn(ji,jj,jk,jn)                                   
317                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
318                  !
319                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn
320                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
321               END DO
322           END DO
323         END DO
324         !
325      END DO
326      !
327   END SUBROUTINE tra_nxt_fix
328
329
330   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )
331      !!----------------------------------------------------------------------
332      !!                   ***  ROUTINE tra_nxt_vvl  ***
333      !!
334      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
335      !!                and swap the tracer fields.
336      !!
337      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
338      !!              - swap tracer fields to prepare the next time_step.
339      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
340      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
341      !!             tn  = ta
342      !!
343      !! ** Action  : - tsb & tsn ready for the next time step
344      !!----------------------------------------------------------------------
345      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
346      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
347      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
348      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
349      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
350      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
351      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
352      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
353      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
354      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
355      !
356      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
357      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
358      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
359      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d, zscale  !   -      -
360      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
361      !!----------------------------------------------------------------------
362      !
363      IF( kt == kit000 )  THEN
364         IF(lwp) WRITE(numout,*)
365         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
366         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
367      ENDIF
368      !
369      IF( cdtype == 'TRA' )  THEN   
370         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
371         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
372         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
373      ELSE                          ! passive tracers case
374         ll_traqsr  = .FALSE.          ! NO solar penetration
375         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
376         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
377      ENDIF
378      !
379      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
380         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
381         ztrd_atf(:,:,:,:) = 0.0_wp
382      ENDIF
383      zfact = 1._wp / p2dt
384      zfact1 = atfp * p2dt
385      zfact2 = zfact1 * r1_rau0
386      DO jn = 1, kjpt     
387         DO jk = 1, jpkm1
388            DO jj = 2, jpjm1
389               DO ji = fs_2, fs_jpim1
390                  ze3t_b = e3t_b(ji,jj,jk)
391                  ze3t_n = e3t_n(ji,jj,jk)
392                  ze3t_a = e3t_a(ji,jj,jk)
393                  !                                         ! tracer content at Before, now and after
394                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
395                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
396                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
397                  !
398                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
399                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
400                  !
401                  ze3t_f = ze3t_n + atfp * ze3t_d
402                  ztc_f  = ztc_n  + atfp * ztc_d
403                  !
404                  zscale = zfact2 * e3t_n(ji,jj,jk) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
405                  ze3t_f = ze3t_f - zscale * ( emp_b(ji,jj) - emp(ji,jj) )
406                  IF ( ll_rnf ) ze3t_f = ze3t_f + zscale * (    rnf_b(ji,jj) -    rnf(ji,jj) )
407                  IF ( ll_isf ) ze3t_f = ze3t_f - zscale * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) )
408
409                  IF( jk == mikt(ji,jj) ) THEN           ! first level
410                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
411                  ENDIF
412                  !
413                  ! solar penetration (temperature only)
414                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
415                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
416                     !
417                  ! river runoff
418                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
419                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
420                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
421                     !
422                  ! ice shelf
423                  IF( ll_isf ) THEN
424                     ! level fully include in the Losch_2008 ice shelf boundary layer
425                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
426                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
427                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
428                     ! level partially include in Losch_2008 ice shelf boundary layer
429                     IF ( jk == misfkb(ji,jj) )                                                   &
430                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
431                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
432                  END IF
433                  !
434                  ze3t_f = 1.e0 / ze3t_f
435                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
436                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
437                  !
438                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
439                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
440                  ENDIF
441                  !
442               END DO
443            END DO
444         END DO
445         !
446      END DO
447      !
448      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
449         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN
450            CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
451            CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
452         ENDIF
453         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
454            DO jn = 1, kjpt
455               CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
456            END DO
457         ENDIF
458         DEALLOCATE( ztrd_atf )
459      ENDIF
460      !
461   END SUBROUTINE tra_nxt_vvl
462
463   !!======================================================================
464END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.