Changeset 4680
- Timestamp:
- 2014-06-20T12:46:45+02:00 (10 years ago)
- Location:
- branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90
r3625 r4680 86 86 zdiv0(:, 1 ) = 0._wp 87 87 zdiv0(:,jpj) = 0._wp 88 IF( .NOT.lk_vopt_loop ) THEN 89 zflu (jpi,:) = 0._wp 90 zflv (jpi,:) = 0._wp 91 zdiv0(1, :) = 0._wp 92 zdiv0(jpi,:) = 0._wp 93 ENDIF 88 zflu (jpi,:) = 0._wp 89 zflv (jpi,:) = 0._wp 90 zdiv0(1, :) = 0._wp 91 zdiv0(jpi,:) = 0._wp 94 92 95 93 zconv = 1._wp !== horizontal diffusion using a Crant-Nicholson scheme ==! -
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90
r4147 r4680 14 14 !! 'key_lim2' : LIM 2.0 sea-ice model 15 15 !!---------------------------------------------------------------------- 16 !!----------------------------------------------------------------------17 16 !! lim_istate_2 : Initialisation of diagnostics ice variables 18 17 !! lim_istate_init_2 : initialization of ice state and namelist read … … 34 33 PUBLIC lim_istate_2 ! routine called by lim_init_2.F90 35 34 36 ! !! ** namelist (namiceini) **37 LOGICAL :: ln_limini ! :Ice initialization state35 ! !! ** namelist (namiceini) ** 36 LOGICAL :: ln_limini ! Ice initialization state 38 37 REAL(wp) :: ttest ! threshold water temperature for initial sea ice 39 38 REAL(wp) :: hninn ! initial snow thickness in the north … … 51 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 52 51 !!---------------------------------------------------------------------- 53 54 52 CONTAINS 55 53 … … 71 69 IF( .NOT. ln_limini ) THEN 72 70 73 tfu(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius]71 tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 74 72 75 73 DO jj = 1, jpj -
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r4333 r4680 83 83 zdiv0(:, 1 ) = 0._wp 84 84 zdiv0(:,jpj) = 0._wp 85 IF( .NOT.lk_vopt_loop ) THEN 86 zflu (jpi,:) = 0._wp 87 zflv (jpi,:) = 0._wp 88 zdiv0(1, :) = 0._wp 89 zdiv0(jpi,:) = 0._wp 90 ENDIF 85 zflu (jpi,:) = 0._wp 86 zflv (jpi,:) = 0._wp 87 zdiv0(1, :) = 0._wp 88 zdiv0(jpi,:) = 0._wp 91 89 92 90 zconv = 1._wp !== horizontal diffusion using a Crant-Nicholson scheme ==! -
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4335 r4680 36 36 PUBLIC lim_istate ! routine called by lim_init.F90 37 37 38 !! * Module variables39 38 ! !!** init namelist (namiceini) ** 40 39 REAL(wp) :: ttest ! threshold water temperature for initial sea ice … … 53 52 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 54 53 !!---------------------------------------------------------------------- 55 56 54 CONTAINS 57 55 … … 121 119 122 120 ! Basal temperature is set to the freezing point of seawater in Celsius 123 t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius]121 t_bo(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 124 122 125 123 DO jj = 1, jpj ! ice if sst <= t-freez + ttest -
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r4619 r4680 551 551 DO jj = 1, jpjm1 552 552 DO ji = 1, jpim1 553 mgrhu(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )554 mgrhv(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )553 mgrhu(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) 554 mgrhv(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) 555 555 END DO 556 556 END DO … … 558 558 DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point 559 559 DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0) 560 e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj )), fse3u_0(ji,jj,mbkt(ji,jj)) )561 e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )560 e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj )), e3u_0(ji,jj,mbkt(ji,jj)) ) 561 e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) 562 562 END DO 563 563 END DO -
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r4619 r4680 2 2 !!============================================================================== 3 3 !! *** MODULE tranpc *** 4 !! Ocean active tracers: non penetrative convecti onscheme4 !! Ocean active tracers: non penetrative convective adjustment scheme 5 5 !!============================================================================== 6 6 !! History : 1.0 ! 1990-09 (G. Madec) Original code … … 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 12 !!---------------------------------------------------------------------- 12 13 … … 14 15 !! tra_npc : apply the non penetrative convection scheme 15 16 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and active tracers 17 USE oce ! ocean dynamics and active tracers 17 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 18 20 USE zdf_oce ! ocean vertical physics 19 21 USE trd_oce ! ocean active tracer trends 20 22 USE trdtra ! ocean active tracer trends 21 USE eosbn2 ! equation of state (eos routine) 23 USE eosbn2 ! equation of state (eos routine) 22 24 USE lbclnk ! lateral boundary conditions (or mpp link) 23 25 USE in_out_manager ! I/O manager … … 29 31 PRIVATE 30 32 31 PUBLIC tra_npc 33 PUBLIC tra_npc ! routine called by step.F90 32 34 33 35 !! * Substitutions 34 36 # include "domzgr_substitute.h90" 35 !!---------------------------------------------------------------------- 36 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 37 !! $Id$ 37 # include "vectopt_loop_substitute.h90" 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 40 !! $Id$ 38 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 39 42 !!---------------------------------------------------------------------- … … 44 47 !! *** ROUTINE tranpc *** 45 48 !! 46 !! ** Purpose : Non penetrative convective adjustment scheme. solve49 !! ** Purpose : Non-penetrative convective adjustment scheme. solve 47 50 !! the static instability of the water column on after fields 48 51 !! while conserving heat and salt contents. 49 52 !! 50 !! ** Method : The algorithm used converges in a maximium of jpk 51 !! iterations. instabilities are treated when the vertical density 52 !! gradient is less than 1.e-5. 53 !! l_trdtra=T: the trend associated with this algorithm is saved. 53 !! ** Method : updated algorithm able to deal with non-linear equation of state 54 !! (i.e. static stability computed locally) 54 55 !! 55 56 !! ** Action : - (ta,sa) after the application od the npc scheme 56 !! - s ave the associated trends (ttrd,strd) ('key_trdtra')57 !! - send the associated trends for on-line diagnostics (l_trdtra=T) 57 58 !! 58 !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.59 !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. 59 60 !!---------------------------------------------------------------------- 60 !61 61 INTEGER, INTENT(in) :: kt ! ocean time-step index 62 62 ! 63 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 64 INTEGER :: inpcc ! number of statically instable water column 65 INTEGER :: inpci ! number of iteration for npc scheme 66 INTEGER :: jiter, jkdown, jkp ! ??? 67 INTEGER :: ikbot, ik, ikup, ikdown ! ??? 68 REAL(wp) :: ze3tot, zta, zsa, zraua, ze3dwn 69 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx, zwy, zwz 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds, zrhop 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. 71 81 !!---------------------------------------------------------------------- 72 82 ! 73 83 IF( nn_timing == 1 ) CALL timing_start('tra_npc') 74 84 ! 75 CALL wrk_alloc(jpi, jpj, jpk, zrhop )76 CALL wrk_alloc(jpi, jpk, zwx, zwy, zwz )77 !78 85 IF( MOD( kt, nn_npc ) == 0 ) THEN 79 80 inpcc = 081 inpci = 082 83 CALL eos( tsa, rhd, zrhop, fsdept_n(:,:,:) ) ! Potential density84 85 IF( l_trdtra ) THEN !* Save ta and sa trends86 ! 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 86 93 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 87 94 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 89 96 ENDIF 90 97 91 ! ! =============== 92 DO jj = 1, jpj ! Vertical slab 93 ! ! =============== 94 ! Static instability pointer 95 ! ---------------------------- 96 DO jk = 1, jpkm1 97 DO ji = 1, jpi 98 zwx(ji,jk) = ( zrhop(ji,jj,jk) - zrhop(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 99 END DO 100 END DO 101 102 ! 1.1 do not consider the boundary points 103 104 ! even if east-west cyclic b. c. do not considere ji=1 or jpi 105 DO jk = 1, jpkm1 106 zwx( 1 ,jk) = 0.e0 107 zwx(jpi,jk) = 0.e0 108 END DO 109 ! even if south-symmetric b. c. used, do not considere jj=1 110 IF( jj == 1 ) zwx(:,:) = 0.e0 111 112 DO jk = 1, jpkm1 113 DO ji = 1, jpi 114 zwx(ji,jk) = 1. 115 IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0 116 END DO 117 END DO 118 119 zwy(:,1) = 0.e0 120 DO ji = 1, jpi 121 DO jk = 1, jpkm1 122 zwy(ji,1) = zwy(ji,1) + zwx(ji,jk) 123 END DO 124 END DO 125 126 zwz(1,1) = 0.e0 127 DO ji = 1, jpi 128 zwz(1,1) = zwz(1,1) + zwy(ji,1) 129 END DO 130 131 inpcc = inpcc + NINT( zwz(1,1) ) 132 133 134 ! 2. Vertical mixing for each instable portion of the density profil 135 ! ------------------------------------------------------------------ 136 137 IF( zwz(1,1) /= 0.e0 ) THEN ! -->> the density profil is statically instable : 138 DO ji = 1, jpi 139 IF( zwy(ji,1) /= 0.e0 ) THEN 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 ) 140 155 ! 141 ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level 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 142 289 ! 143 DO jiter = 1, jpk ! vertical iteration 144 ! 145 ! search of ikup : the first static instability from the sea surface 146 ! 147 ik = 0 148 220 CONTINUE 149 ik = ik + 1 150 IF( ik >= ikbot ) GO TO 200 151 zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1) 152 IF( zwx(ji,ik) <= 0.e0 ) GO TO 220 153 ikup = ik 154 ! the density profil is instable below ikup 155 ! ikdown : bottom of the instable portion of the density profil 156 ! search of ikdown and vertical mixing from ikup to ikdown 157 ! 158 ze3tot= fse3t(ji,jj,ikup) 159 zta = tsa (ji,jj,ikup,jp_tem) 160 zsa = tsa (ji,jj,ikup,jp_sal) 161 zraua = zrhop(ji,jj,ikup) 162 ! 163 DO jkdown = ikup+1, ikbot-1 164 IF( zraua <= zrhop(ji,jj,jkdown) ) THEN 165 ikdown = jkdown 166 GO TO 240 167 ENDIF 168 ze3dwn = fse3t(ji,jj,jkdown) 169 ze3tot = ze3tot + ze3dwn 170 zta = ( zta*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_tem)*ze3dwn )/ze3tot 171 zsa = ( zsa*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_sal)*ze3dwn )/ze3tot 172 zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot 173 inpci = inpci+1 174 END DO 175 ikdown = ikbot-1 176 240 CONTINUE 177 ! 178 DO jkp = ikup, ikdown-1 179 tsa (ji,jj,jkp,jp_tem) = zta 180 tsa (ji,jj,jkp,jp_sal) = zsa 181 zrhop(ji,jj,jkp ) = zraua 182 END DO 183 IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN 184 tsa (ji,jj,jkp,jp_tem) = zta 185 tsa (ji,jj,jkp,jp_sal) = zsa 186 zrhop(ji,jj,ikdown ) = zraua 187 ENDIF 188 END DO 189 ENDIF 190 200 CONTINUE 191 END DO 192 ! <<-- no more static instability on slab jj 193 ENDIF 194 ! ! =============== 195 END DO ! End of slab 196 ! ! =============== 197 ! 198 IF( l_trdtra ) THEN ! save the Non penetrative mixing trends for diagnostic 199 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 200 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 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 201 315 CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 202 316 CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 203 317 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 204 318 ENDIF 205 206 ! Lateral boundary conditions on ( ta, sa ) ( Unchanged sign) 207 ! ------------------------------============ 319 ! 208 320 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 209 210 211 ! 2. non penetrative convective scheme statistics212 ! -----------------------------------------------213 IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN214 IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable', &215 & ' water column : ',inpcc, ' number of iteration : ',inpci216 ENDIF217 !218 ENDIF219 !220 CALL wrk_dealloc(jpi, jpj, jpk, zrhop )221 CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz )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 222 334 ! 223 335 IF( nn_timing == 1 ) CALL timing_stop('tra_npc')
Note: See TracChangeset
for help on using the changeset viewer.