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.
tranpc.F90 in NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/TRA/tranpc.F90 @ 12165

Last change on this file since 12165 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convective adjustment scheme
5   !!==============================================================================
6   !! History :  1.0  ! 1990-09  (G. Madec)  Original code
7   !!                 ! 1996-01  (G. Madec)  statement function for e3
8   !!   NEMO     1.0  ! 2002-06  (G. Madec)  free form F90
9   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
10   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
11   !!            3.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_npc       : apply the non penetrative convection scheme
16   !!----------------------------------------------------------------------
17   USE oce            ! ocean dynamics and active tracers
18   USE dom_oce        ! ocean space and time domain
19   USE phycst         ! physical constants
20   USE zdf_oce        ! ocean vertical physics
21   USE trd_oce        ! ocean active tracer trends
22   USE trdtra         ! ocean active tracer trends
23   USE eosbn2         ! equation of state (eos routine)
24   !
25   USE lbclnk         ! lateral boundary conditions (or mpp link)
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_npc    ! routine called by step.F90
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
39   !! $Id$
40   !! Software governed by the CeCILL license (see ./LICENSE)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE tra_npc( kt )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE tranpc  ***
47      !!
48      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
49      !!      the static instability of the water column on after fields
50      !!      while conserving heat and salt contents.
51      !!
52      !! ** Method  : updated algorithm able to deal with non-linear equation of state
53      !!              (i.e. static stability computed locally)
54      !!
55      !! ** Action  : - tsa: after tracers with the application of the npc scheme
56      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
57      !!
58      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
59      !!----------------------------------------------------------------------
60      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
61      !
62      INTEGER  ::   ji, jj, jk   ! dummy loop indices
63      INTEGER  ::   inpcc        ! number of statically instable water column
64      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers
65      LOGICAL  ::   l_bottom_reached, l_column_treated
66      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
67      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
68      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0)
69      REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point...
70      REAL(wp), DIMENSION(        jpk,jpts) ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
71      REAL(wp), DIMENSION(jpi,jpj,jpk     ) ::   zn2          ! N^2
72      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zab          ! alpha and beta
73      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
74      !
75      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
76      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
77      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
78      !!----------------------------------------------------------------------
79      !
80      IF( ln_timing )   CALL timing_start('tra_npc')
81      !
82      IF( MOD( kt, nn_npc ) == 0 ) THEN
83         !
84         IF( l_trdtra )   THEN                    !* Save initial after fields
85            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
86            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
87            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
88         ENDIF
89         !
90         IF( l_LB_debug ) THEN
91            ! Location of 1 known convection site to follow what's happening in the water column
92            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
93            nncpu = 1  ;            ! the CPU domain contains the convection spot
94            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
95         ENDIF
96         !
97         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points)
98         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points)
99         !
100         inpcc = 0
101         !
102         DO jj = 2, jpjm1                 ! interior column only
103            DO ji = fs_2, fs_jpim1
104               !
105               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
106                  !                                     ! consider one ocean column
107                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
108                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
109                  !
110                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
111                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
112                  zvn2(:)         = zn2(ji,jj,:)            ! N^2
113                  !
114                  IF( l_LB_debug ) THEN                  !LB debug:
115                     lp_monitor_point = .FALSE.
116                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
117                     ! writing only if on CPU domain where conv region is:
118                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                     
119                  ENDIF                                  !LB debug  end
120                  !
121                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
122                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2)
123                  ilayer = 0
124                  jiter  = 0
125                  l_column_treated = .FALSE.
126                  !
127                  DO WHILE ( .NOT. l_column_treated )
128                     !
129                     jiter = jiter + 1
130                     !
131                     IF( jiter >= 400 ) EXIT
132                     !
133                     l_bottom_reached = .FALSE.
134                     !
135                     DO WHILE ( .NOT. l_bottom_reached )
136                        !
137                        ikp = ikp + 1
138                        !
139                        !! Testing level ikp for instability
140                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
142                           !
143                           ilayer = ilayer + 1    ! yet another instable portion of the water column found....
144                           !
145                           IF( lp_monitor_point ) THEN
146                              WRITE(numout,*)
147                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
148                                 WRITE(numout,*)
149                                 WRITE(numout,*) 'Time step = ',kt,' !!!'
150                              ENDIF
151                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   &
152                                 &                                    ' in column! Starting at ikp =', ikp
153                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj
154                              DO jk = 1, klc1
155                                 WRITE(numout,*) jk, zvn2(jk)
156                              END DO
157                              WRITE(numout,*)
158                           ENDIF
159                           !
160                           IF( jiter == 1 )   inpcc = inpcc + 1 
161                           !
162                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
163                           !
164                           !! ikup is the uppermost point where mixing will start:
165                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
166                           !
167                           !! If the points above ikp-1 have N2 == 0 they must also be mixed:
168                           IF( ikp > 2 ) THEN
169                              DO jk = ikp-1, 2, -1
170                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN
171                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
172                                 ELSE
173                                    EXIT
174                                 ENDIF
175                              END DO
176                           ENDIF
177                           !
178                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
179                           !
180                           zsum_temp = 0._wp
181                           zsum_sali = 0._wp
182                           zsum_alfa = 0._wp
183                           zsum_beta = 0._wp
184                           zsum_z    = 0._wp
185                                                   
186                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
187                              !
188                              zdz       = e3t_n(ji,jj,jk)
189                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
190                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
191                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
192                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
193                              zsum_z    = zsum_z    + zdz
194                              !                             
195                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
196                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
197                              IF( zvn2(jk+1) > zn2_zero ) EXIT
198                           END DO
199                         
200                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
201                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
202
203                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
204                           zta   = zsum_temp/zsum_z
205                           zsa   = zsum_sali/zsum_z
206                           zalfa = zsum_alfa/zsum_z
207                           zbeta = zsum_beta/zsum_z
208
209                           IF( lp_monitor_point ) THEN
210                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
211                                 &            ' and ikdown =',ikdown,', in layer #',ilayer
212                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
213                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
214                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
215                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
216                           ENDIF
217
218                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
219                           DO jk = ikup, ikdown
220                              zvts(jk,jp_tem) = zta
221                              zvts(jk,jp_sal) = zsa
222                              zvab(jk,jp_tem) = zalfa
223                              zvab(jk,jp_sal) = zbeta
224                           END DO
225                           
226                           
227                           !! Updating N2 in the relvant portion of the water column
228                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
229                           !! => Need to re-compute N2! will use Alpha and Beta!
230                           
231                           ikup   = MAX(2,ikup)         ! ikup can never be 1 !
232                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
233                           
234                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
235
236                              !! Interpolating alfa and beta at W point:
237                              zrw =  (gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk)) &
238                                 & / (gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk))
239                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
240                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
241
242                              !! N2 at W point, doing exactly as in eosbn2.F90:
243                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
244                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
245                                 &       / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
246
247                              !! OR, faster  => just considering the vertical gradient of density
248                              !! as only the signa maters...
249                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
250                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
251
252                           END DO
253                       
254                           ikp = MIN(ikdown+1,ikbot)
255                           
256
257                        ENDIF  !IF( zvn2(ikp) < 0. )
258
259
260                        IF( ikp == ikbot ) l_bottom_reached = .TRUE.
261                        !
262                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
263
264                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
265                   
266                     ! ******* At this stage ikp == ikbot ! *******
267                   
268                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
269                        !
270                        IF( lp_monitor_point ) THEN
271                           WRITE(numout,*)
272                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
273                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
274                           DO jk = 1, klc1
275                              WRITE(numout,*) jk, zvn2(jk)
276                           END DO
277                           WRITE(numout,*)
278                        ENDIF
279                        !
280                        ikp    = 1     ! starting again at the surface for the next iteration
281                        ilayer = 0
282                     ENDIF
283                     !
284                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
285                     !
286                  END DO ! DO WHILE ( .NOT. l_column_treated )
287
288                  !! Updating tsa:
289                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
290                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
291
292                  !! LB:  Potentially some other global variable beside theta and S can be treated here
293                  !!      like BGC tracers.
294
295                  IF( lp_monitor_point )   WRITE(numout,*)
296
297               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
298
299            END DO ! ji
300         END DO ! jj
301         !
302         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
303            z1_r2dt = 1._wp / (2._wp * rdt)
304            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
305            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
306            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
307            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
308            DEALLOCATE( ztrdt, ztrds )
309         ENDIF
310         !
311         CALL lbc_lnk_multi( 'tranpc', tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1. )
312         !
313         IF( lwp .AND. l_LB_debug ) THEN
314            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
315            WRITE(numout,*)
316         ENDIF
317         !
318      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
319      !
320      IF( ln_timing )   CALL timing_stop('tra_npc')
321      !
322   END SUBROUTINE tra_npc
323
324   !!======================================================================
325END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.