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.
Changeset 5386 – NEMO

Changeset 5386


Ignore:
Timestamp:
2015-06-09T16:05:14+02:00 (9 years ago)
Author:
gm
Message:

#1092 minor bug correction on tranpc.F90 algorithm

Location:
trunk/NEMOGCM/NEMO/OPA_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r4990 r5386  
    99   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90 
    1010   !!            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. 
     11   !!            3.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 
    1212   !!---------------------------------------------------------------------- 
    1313 
     
    6464      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6565      INTEGER  ::   inpcc        ! number of statically instable water column 
    66       INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers 
     66      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers 
    6767      LOGICAL  ::   l_bottom_reached, l_column_treated 
    6868      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
    6969      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     70      REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0) 
    7071      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
    7172      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
     
    7576      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
    7677      ! 
    77       !!LB debug: 
    78       LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 
    79       INTEGER :: ilc1, jlc1, klc1, nncpu 
    80       LOGICAL :: lp_monitor_point = .FALSE. 
    81       !!LB debug. 
     78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 
     79      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1" 
     80      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu" 
    8281      !!---------------------------------------------------------------------- 
    8382      ! 
     
    9796         ENDIF 
    9897 
    99          !LB debug: 
    100          IF( lwp .AND. l_LB_debug ) THEN 
    101             WRITE(numout,*) 
    102             WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea 
    103          ENDIF 
    104          !LBdebug: Monitoring of 1 column subject to convection... 
    10598         IF( l_LB_debug ) THEN 
    106             ! Location of 1 known convection spot to follow what's happening in the water column 
    107             ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 
    108             nncpu = 15  ; ! the CPU domain contains the convection spot 
    109             !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 
    110             !nncpu = 54  ; ! the CPU domain contains the convection spot 
     99            ! Location of 1 known convection site to follow what's happening in the water column 
     100            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...            
     101            nncpu = 1  ;            ! the CPU domain contains the convection spot 
    111102            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
    112103         ENDIF 
    113          !LBdebug. 
    114  
    115          CALL eos_rab( tsa, zab )         ! after alpha and beta 
    116          CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala 
     104          
     105         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points) 
     106         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points) 
    117107         
    118108         inpcc = 0 
     
    134124                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
    135125                     ! writing only if on CPU domain where conv region is: 
    136                      lp_monitor_point = (narea == nncpu).AND.lp_monitor_point  
    137                       
    138                      IF(lp_monitor_point) THEN 
    139                         WRITE(numout,*) '' ;WRITE(numout,*) '' ; 
    140                        WRITE(numout,'("Time step = ",i6.6," !!!")')  kt 
    141                         WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj 
    142                         DO jk = 1, klc1 
    143                            WRITE(numout,*) jk, zvn2(jk) 
    144                         END DO 
    145                         WRITE(numout,*) ' ' 
    146                      ENDIF 
     126                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                       
    147127                  ENDIF                                  !LB debug  end 
    148128 
    149129                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
    150                   ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2) 
     130                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2) 
    151131                  ilayer = 0 
    152132                  jiter  = 0 
     
    163143                     DO WHILE ( .NOT. l_bottom_reached ) 
    164144 
    165                         ik = ik + 1 
     145                        ikp = ikp + 1 
    166146                        
    167                         !! Checking level ik for instability 
     147                        !! Testing level ikp for instability 
    168148                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    169  
    170                         IF( zvn2(ik) < 0. ) THEN ! Instability found! 
    171  
    172                            ikm  = ik              ! first level whith negative N2 
    173                            ilayer = ilayer + 1    ! yet another layer found.... 
    174                            IF(jiter == 1) inpcc = inpcc + 1 
    175  
    176                            IF(l_LB_debug .AND. lp_monitor_point) & 
    177                               & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 
    178                               & ' inpcc =', inpcc 
    179  
    180                            !! Case we mix with upper regions where N2==0: 
    181                            !! All the points above ikup where N2 == 0 must also be mixed => we go 
    182                            !! upward to find a new ikup, where the layer doesn't have N2==0 
    183                            ikup = ikm 
    184                            DO jk = ikm, 2, -1 
    185                               ikup = ikup - 1 
    186                               IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 
    187                            END DO 
    188                            
    189                            ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 
    190                            IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 
    191  
    192                            
    193                            IF( lp_monitor_point )   WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 
    194                            
     149                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found! 
     150 
     151                           ilayer = ilayer + 1    ! yet another instable portion of the water column found.... 
     152 
     153                           IF( lp_monitor_point ) THEN  
     154                              WRITE(numout,*) 
     155                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability 
     156                                 WRITE(numout,*) 
     157                                 WRITE(numout,*) 'Time step = ',kt,' !!!' 
     158                              ENDIF 
     159                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   & 
     160                                 &                                    ' in column! Starting at ikp =', ikp 
     161                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj 
     162                              DO jk = 1, klc1 
     163                                 WRITE(numout,*) jk, zvn2(jk) 
     164                              END DO 
     165                              WRITE(numout,*) 
     166                           ENDIF 
     167                            
     168 
     169                           IF( jiter == 1 )   inpcc = inpcc + 1  
     170 
     171                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 
     172 
     173                           !! ikup is the uppermost point where mixing will start: 
     174                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying 
     175                            
     176                           !! If the points above ikp-1 have N2 == 0 they must also be mixed: 
     177                           IF( ikp > 2 ) THEN 
     178                              DO jk = ikp-1, 2, -1 
     179                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN 
     180                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing 
     181                                 ELSE 
     182                                    EXIT 
     183                                 ENDIF 
     184                              END DO 
     185                           ENDIF 
     186                            
     187                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1') 
     188 
    195189                           zsum_temp = 0._wp 
    196190                           zsum_sali = 0._wp 
     
    199193                           zsum_z    = 0._wp 
    200194                                                     
    201                            DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column 
    202                               ! 
    203                               IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk 
     195                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column 
    204196                              ! 
    205197                              zdz       = fse3t(ji,jj,jk) 
     
    209201                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
    210202                              zsum_z    = zsum_z    + zdz 
    211                               ! 
    212                               !! EXIT if we found the bottom of the unstable portion of the water column     
    213                               IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT 
     203                              !                               
     204                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 
     205                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 
     206                              IF( zvn2(jk+1) > zn2_zero ) EXIT 
    214207                           END DO 
    215208                           
    216                            !ik     = jk !LB remove? 
    217                            ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 
    218                            
    219                            IF(l_LB_debug .AND. lp_monitor_point) & 
    220                               &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer 
    221                            
    222                            ! Mixing Temperature and salinity between ikup and ikdown: 
     209                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 
     210                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2') 
     211 
     212                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included: 
    223213                           zta   = zsum_temp/zsum_z 
    224214                           zsa   = zsum_sali/zsum_z 
     
    226216                           zbeta = zsum_beta/zsum_z 
    227217 
    228                            IF(l_LB_debug .AND. lp_monitor_point) THEN 
     218                           IF( lp_monitor_point ) THEN 
     219                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   & 
     220                                 &            ' and ikdown =',ikdown,', in layer #',ilayer 
    229221                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
    230222                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
    231                               WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa 
     223                              WRITE(numout,*) '  => Mean Alfa in that portion =', zalfa 
    232224                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
    233225                           ENDIF 
     
    240232                              zvab(jk,jp_sal) = zbeta 
    241233                           END DO 
    242                            ! 
    243                            !! Before updating N2, it is possible that another unstable 
    244                            !! layer exists underneath the one we just homogeneized! 
    245                            ik = ikdown 
    246                            !  
    247                         ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN 
    248                         ! 
    249                         IF( ik == ikbot ) l_bottom_reached = .TRUE. 
     234                            
     235                            
     236                           !! Updating N2 in the relvant portion of the water column 
     237                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     238                           !! => Need to re-compute N2! will use Alpha and Beta! 
     239                            
     240                           ikup   = MAX(2,ikup)         ! ikup can never be 1 ! 
     241                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 
     242                            
     243                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown! 
     244 
     245                              !! Interpolating alfa and beta at W point: 
     246                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     247                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     248                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     249                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     250 
     251                              !! N2 at W point, doing exactly as in eosbn2.F90: 
     252                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     253                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     254                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     255 
     256                              !! OR, faster  => just considering the vertical gradient of density 
     257                              !! as only the signa maters... 
     258                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     259                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  ) 
     260 
     261                           END DO 
     262                         
     263                           ikp = MIN(ikdown+1,ikbot) 
     264                            
     265 
     266                        ENDIF  !IF( zvn2(ikp) < 0. ) 
     267 
     268 
     269                        IF( ikp == ikbot ) l_bottom_reached = .TRUE. 
    250270                        ! 
    251271                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
    252272 
    253                      IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1' 
     273                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3') 
    254274                     
    255                      ! ******* At this stage ik == ikbot ! ******* 
     275                     ! ******* At this stage ikp == ikbot ! ******* 
    256276                     
    257                      IF( ilayer > 0 ) THEN 
    258                         !! least an unstable layer has been found 
    259                         !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
    260                         !! => Need to re-compute N2! will use Alpha and Beta! 
     277                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found 
    261278                        ! 
    262                         DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!      
    263                            !! Doing exactly as in eosbn2.F90: 
    264                            !! * Except that we only are interested in the sign of N2 !!! 
    265                            !!   => just considering the vertical gradient of density 
    266                            zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
    267                                & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
    268                            zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
    269                            zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
    270                            
    271                            !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
    272                            !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
    273                            !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
    274                            zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
    275                                 &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )   
    276                         END DO 
    277  
    278                         IF(l_LB_debug .AND. lp_monitor_point) THEN 
    279                            WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 
    280                               & jiter, ji,jj 
     279                        IF( lp_monitor_point ) THEN 
     280                           WRITE(numout,*) 
     281                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)' 
     282                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:' 
    281283                           DO jk = 1, klc1 
    282284                              WRITE(numout,*) jk, zvn2(jk) 
    283285                           END DO 
    284                            WRITE(numout,*) ' ' 
     286                           WRITE(numout,*) 
    285287                        ENDIF 
    286  
    287                         ik     = 1  ! starting again at the surface for the next iteration 
     288                        ! 
     289                        ikp    = 1     ! starting again at the surface for the next iteration 
    288290                        ilayer = 0 
    289291                     ENDIF 
    290292                     ! 
    291                      IF( ik >= ikbot ) THEN 
    292                         IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---' 
    293                         l_column_treated = .TRUE. 
    294                      ENDIF 
     293                     IF( ikp >= ikbot )   l_column_treated = .TRUE. 
    295294                     ! 
    296295                  END DO ! DO WHILE ( .NOT. l_column_treated ) 
     
    300299                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 
    301300 
    302                   !! lolo:  Should we update something else???? 
    303                   !! => like alpha and beta? 
    304  
    305                   IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '' 
     301                  !! LB:  Potentially some other global variable beside theta and S can be treated here 
     302                  !!      like BGC tracers. 
     303 
     304                  IF( lp_monitor_point )   WRITE(numout,*) 
    306305 
    307306               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN 
     
    321320         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    322321         ! 
    323          IF(lwp) THEN 
    324             WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 
    325             WRITE(numout,*)' => number of statically instable water column : ',inpcc 
    326             WRITE(numout,*) '' ; WRITE(numout,*) '' 
     322         IF( lwp .AND. l_LB_debug ) THEN 
     323            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc 
     324            WRITE(numout,*) 
    327325         ENDIF 
    328326         ! 
  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r5120 r5386  
    124124      IF(lwp) WRITE(numout,*) '   convection :' 
    125125      ! 
    126       IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: non penetrative convective scheme is not working',   & 
    127          &                                       ' set ln_zdfnpc to FALSE' ) 
     126#if defined key_top 
     127      IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_init: npc scheme is not working with key_top' ) 
     128#endif 
    128129      ! 
    129130      ioptio = 0 
Note: See TracChangeset for help on using the changeset viewer.