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 branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 4680

Last change on this file since 4680 was 4680, checked in by gm, 10 years ago

#1350 : new version of tranpc.F90 based local static stability
  • Property svn:keywords set to Id
File size: 16.9 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.7  ! 2014-06  (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   USE lbclnk          ! lateral boundary conditions (or mpp link)
25   USE in_out_manager  ! I/O manager
26   USE lib_mpp         ! MPP library
27   USE wrk_nemo        ! Memory Allocation
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 "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tra_npc( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tranpc  ***
48      !!
49      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
50      !!      the static instability of the water column on after fields
51      !!      while conserving heat and salt contents.
52      !!
53      !! ** Method  : updated algorithm able to deal with non-linear equation of state
54      !!              (i.e. static stability computed locally)
55      !!
56      !! ** Action  : - (ta,sa) after the application od the npc scheme
57      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
58      !!
59      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
60      !!----------------------------------------------------------------------
61      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
62      !
63      INTEGER  ::   ji, jj, jk   ! dummy loop indices
64      INTEGER  ::   inpcc        ! number of statically instable water column
65      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers
66      LOGICAL  ::   l_bottom_reached, l_column_treated
67      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
68      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
69      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point...
70      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point...
71      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta
72      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2
73      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta
74      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace
75      !
76      !!LB debug:
77      LOGICAL, PARAMETER :: l_LB_debug = .FALSE.
78      INTEGER :: ilc1, jlc1, klc1, nncpu
79      LOGICAL :: lp_monitor_point
80      !!LB debug.
81      !!----------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('tra_npc')
84      !
85      IF( MOD( kt, nn_npc ) == 0 ) THEN
86         !
87         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2
88         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta
89         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj
90         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj
91
92         IF( l_trdtra )   THEN                    !* Save initial after fields
93            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
94            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
95            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
96         ENDIF
97
98         !LB debug:
99         IF( lwp .AND. l_LB_debug ) THEN
100            WRITE(numout,*) '' ; WRITE(numout,*) ''
101            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea
102         ENDIF
103         !LBdebug: Monitoring of 1 column subject to convection...
104         IF( l_LB_debug ) THEN
105            ! Location of 1 known convection spot to follow what's happening in the water column
106            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus:
107            nncpu = 15  ; ! the CPU domain contains the convection spot
108            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus:
109            !nncpu = 54  ; ! the CPU domain contains the convection spot
110            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
111         ENDIF
112         !LBdebug.
113
114         CALL eos_rab( tsa, zab )         ! after alpha and beta
115         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala
116       
117         inpcc = 0
118
119         DO jj = 2, jpjm1                 ! interior column only
120            DO ji = fs_2, fs_jpim1
121               !
122               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
123                  !                                     ! consider one ocean column
124                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
125                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
126
127                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
128                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
129                  zvn2(:)         = zn2(ji,jj,:)            ! N^2
130                 
131                  IF( l_LB_debug ) THEN                  !LB debug:
132                     lp_monitor_point = .FALSE.
133                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
134                     ! writing only if on CPU domain where conv region is:
135                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 
136                     
137                     IF(lp_monitor_point) THEN
138                        WRITE(numout,*) '' ;WRITE(numout,*) '' ;
139                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt
140                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') , ji,jj
141                        DO jk = 1, klc1
142                           WRITE(numout,*) jk, zvn2(jk)
143                        END DO
144                        WRITE(numout,*) ' '
145                     ENDIF
146                  ENDIF                                  !LB debug  end
147
148                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
149                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2)
150                  ilayer = 0
151                  jiter  = 0
152                  l_column_treated = .FALSE.
153                 
154                  DO WHILE ( .NOT. l_column_treated )
155                     !
156                     jiter = jiter + 1
157                   
158                     IF( jiter >= 400 ) EXIT
159                   
160                     l_bottom_reached = .FALSE.
161
162                     DO WHILE ( .NOT. l_bottom_reached )
163
164                        ik = ik + 1
165                       
166                        !! Checking level ik for instability
167                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
168
169                        IF( zvn2(ik) < 0. ) THEN ! Instability found!
170
171                           ikm  = ik              ! first level whith negative N2
172                           ilayer = ilayer + 1    ! yet another layer found....
173                           IF(jiter == 1) inpcc = inpcc + 1
174
175                           IF(l_LB_debug .AND. lp_monitor_point) &
176                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, &
177                              & ' inpcc =', inpcc
178
179                           !! Case we mix with upper regions where N2==0:
180                           !! All the points above ikup where N2 == 0 must also be mixed => we go
181                           !! upward to find a new ikup, where the layer doesn't have N2==0
182                           ikup = ikm
183                           DO jk = ikm, 2, -1
184                              ikup = ikup - 1
185                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT
186                           END DO
187                         
188                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0)
189                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ;
190
191                         
192                           IF(lp_monitor_point) WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer
193                         
194                           zsum_temp = 0._wp
195                           zsum_sali = 0._wp
196                           zsum_alfa = 0._wp
197                           zsum_beta = 0._wp
198                           zsum_z    = 0._wp
199                                                   
200                           DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column
201                              !
202                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk
203                              !
204                              zdz       = fse3t(ji,jj,jk)
205                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
206                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
207                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
208                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
209                              zsum_z    = zsum_z    + zdz
210                              !
211                              !! EXIT if we found the bottom of the unstable portion of the water column   
212                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT
213                           END DO
214                         
215                           !ik     = jk !LB remove?
216                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2
217                         
218                           IF(l_LB_debug .AND. lp_monitor_point) &
219                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer
220                         
221                           ! Mixing Temperature and salinity between ikup and ikdown:
222                           zta   = zsum_temp/zsum_z
223                           zsa   = zsum_sali/zsum_z
224                           zalfa = zsum_alfa/zsum_z
225                           zbeta = zsum_beta/zsum_z
226
227                           IF(l_LB_debug .AND. lp_monitor_point) THEN
228                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
229                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
230                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa
231                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
232                           ENDIF
233
234                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
235                           DO jk = ikup, ikdown
236                              zvts(jk,jp_tem) = zta
237                              zvts(jk,jp_sal) = zsa
238                              zvab(jk,jp_tem) = zalfa
239                              zvab(jk,jp_sal) = zbeta
240                           END DO
241                           !
242                           !! Before updating N2, it is possible that another unstable
243                           !! layer exists underneath the one we just homogeneized!
244                           ik = ikdown
245                           !
246                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN
247                        !
248                        IF( ik == ikbot ) l_bottom_reached = .TRUE.
249                        !
250                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
251
252                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1'
253                   
254                     ! ******* At this stage ik == ikbot ! *******
255                   
256                     IF( ilayer > 0 ) THEN
257                        !! least an unstable layer has been found
258                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
259                        !! => Need to re-compute N2! will use Alpha and Beta!
260                        !
261                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!     
262                           !! Doing exactly as in eosbn2.F90:
263                           !! * Except that we only are interested in the sign of N2 !!!
264                           !!   => just considering the vertical gradient of density
265                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
266                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
267                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
268                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
269                         
270                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
271                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
272                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
273                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
274                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  ) 
275                        END DO
276
277                        IF(l_LB_debug .AND. lp_monitor_point) THEN
278                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') &
279                              & jiter, ji,jj
280                           DO jk = 1, klc1
281                              WRITE(numout,*) jk, zvn2(jk)
282                           END DO
283                           WRITE(numout,*) ' '
284                        ENDIF
285
286                        ik     = 1  ! starting again at the surface for the next iteration
287                        ilayer = 0
288                     ENDIF
289                     !
290                     IF( ik >= ikbot ) THEN
291                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---'
292                        l_column_treated = .TRUE.
293                     ENDIF
294                     !
295                  END DO ! DO WHILE ( .NOT. l_column_treated )
296
297                  !! Updating tsa:
298                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
299                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
300
301                  !! lolo:  Should we update something else????
302                  !! => like alpha and beta?
303
304                  IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ''
305
306               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
307
308            END DO ! ji
309         END DO ! jj
310         !
311         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
312            z1_r2dt = 1._wp / (2._wp * rdt)
313            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
314            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
315            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
316            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
317            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
318         ENDIF
319         !
320         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
321         !
322         IF(lwp) THEN
323            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt
324            WRITE(numout,*)' => number of statically instable water column : ',inpcc
325            WRITE(numout,*) '' ; WRITE(numout,*) ''
326         ENDIF
327         !
328         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
329         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
330         CALL wrk_dealloc(jpk, zvn2 )
331         CALL wrk_dealloc(jpk, 2, zvts, zvab )
332         !
333      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
334      !
335      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
336      !
337   END SUBROUTINE tra_npc
338
339   !!======================================================================
340END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.