Changeset 5386
- Timestamp:
- 2015-06-09T16:05:14+02:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 edited
-
TRA/tranpc.F90 (modified) (12 diffs)
-
ZDF/zdfini.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r4990 r5386 9 9 !! 3.0 ! 2008-06 (G. Madec) applied on ta, sa and called before tranxt in step.F90 10 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.11 !! 3.6 ! 2015-05 (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 12 12 !!---------------------------------------------------------------------- 13 13 … … 64 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 65 INTEGER :: inpcc ! number of statically instable water column 66 INTEGER :: jiter, ikbot, ik , ikup, ikdown, ilayer, ikm! local integers66 INTEGER :: jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low ! local integers 67 67 LOGICAL :: l_bottom_reached, l_column_treated 68 68 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 69 69 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) 70 71 REAL(wp), POINTER, DIMENSION(:) :: zvn2 ! vertical profile of N2 at 1 given point... 71 72 REAL(wp), POINTER, DIMENSION(:,:) :: zvts ! vertical profile of T and S at 1 given point... … … 75 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace 76 77 ! 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" 82 81 !!---------------------------------------------------------------------- 83 82 ! … … 97 96 ENDIF 98 97 99 !LB debug:100 IF( lwp .AND. l_LB_debug ) THEN101 WRITE(numout,*)102 WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea103 ENDIF104 !LBdebug: Monitoring of 1 column subject to convection...105 98 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 111 102 klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... 112 103 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) 117 107 118 108 inpcc = 0 … … 134 124 IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 135 125 ! 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 147 127 ENDIF !LB debug end 148 128 149 129 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) 151 131 ilayer = 0 152 132 jiter = 0 … … 163 143 DO WHILE ( .NOT. l_bottom_reached ) 164 144 165 ik = ik+ 1145 ikp = ikp + 1 166 146 167 !! Checking level ikfor instability147 !! Testing level ikp for instability 168 148 !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 195 189 zsum_temp = 0._wp 196 190 zsum_sali = 0._wp … … 199 193 zsum_z = 0._wp 200 194 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 204 196 ! 205 197 zdz = fse3t(ji,jj,jk) … … 209 201 zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 210 202 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 214 207 END DO 215 208 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: 223 213 zta = zsum_temp/zsum_z 224 214 zsa = zsum_sali/zsum_z … … 226 216 zbeta = zsum_beta/zsum_z 227 217 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 229 221 WRITE(numout,*) ' => Mean temp. in that portion =', zta 230 222 WRITE(numout,*) ' => Mean sali. in that portion =', zsa 231 WRITE(numout,*) ' => Mean Al phain that portion =', zalfa223 WRITE(numout,*) ' => Mean Alfa in that portion =', zalfa 232 224 WRITE(numout,*) ' => Mean Beta in that portion =', zbeta 233 225 ENDIF … … 240 232 zvab(jk,jp_sal) = zbeta 241 233 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. 250 270 ! 251 271 END DO ! DO WHILE ( .NOT. l_bottom_reached ) 252 272 253 IF( ik /= ikbot ) STOP 'ERROR: tranpc.F90 => PROBLEM #1'273 IF( ikp /= ikbot ) CALL ctl_stop( 'tra_npc : PROBLEM #3') 254 274 255 ! ******* At this stage ik == ikbot ! *******275 ! ******* At this stage ikp == ikbot ! ******* 256 276 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 261 278 ! 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:' 281 283 DO jk = 1, klc1 282 284 WRITE(numout,*) jk, zvn2(jk) 283 285 END DO 284 WRITE(numout,*) ' '286 WRITE(numout,*) 285 287 ENDIF 286 287 ik = 1! starting again at the surface for the next iteration288 ! 289 ikp = 1 ! starting again at the surface for the next iteration 288 290 ilayer = 0 289 291 ENDIF 290 292 ! 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. 295 294 ! 296 295 END DO ! DO WHILE ( .NOT. l_column_treated ) … … 300 299 tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 301 300 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,*) 306 305 307 306 ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN … … 321 320 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 322 321 ! 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,*) 327 325 ENDIF 328 326 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
r5120 r5386 124 124 IF(lwp) WRITE(numout,*) ' convection :' 125 125 ! 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 128 129 ! 129 130 ioptio = 0
Note: See TracChangeset
for help on using the changeset viewer.
