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/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/tranpc.F90 @ 14088

Last change on this file since 14088 was 14072, checked in by laurent, 4 years ago

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369

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