Changeset 11422 for NEMO/branches/2019/fix_vvl_ticket1791/src/OCE
- Timestamp:
- 2019-08-08T15:40:47+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/fix_vvl_ticket1791/src/OCE
- Files:
-
- 16 edited
-
BDY/bdydta.F90 (modified) (1 diff)
-
BDY/bdydyn2d.F90 (modified) (2 diffs)
-
BDY/bdyice.F90 (modified) (3 diffs)
-
DIA/diawri.F90 (modified) (4 diffs)
-
DOM/dommsk.F90 (modified) (1 diff)
-
DOM/domvvl.F90 (modified) (1 diff)
-
DYN/dynhpg.F90 (modified) (6 diffs)
-
DYN/dynzdf.F90 (modified) (6 diffs)
-
DYN/sshwzv.F90 (modified) (5 diffs)
-
LBC/mppini.F90 (modified) (1 diff)
-
SBC/sbcapr.F90 (modified) (1 diff)
-
STO/stopar.F90 (modified) (1 diff)
-
TRA/traadv_fct.F90 (modified) (8 diffs)
-
TRA/traqsr.F90 (modified) (1 diff)
-
step.F90 (modified) (2 diffs)
-
stpctl.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/BDY/bdydta.F90
r10951 r11422 357 357 jfld_hts = jfld_htst(jbdy) 358 358 jfld_ai = jfld_ait(jbdy) 359 IF ( jpl /= 1 .AND. nice_cat == 1 ) THEN ! case input cat = 1 360 CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 361 & dta_bdy(jbdy)%h_i , dta_bdy(jbdy)%h_s , dta_bdy(jbdy)%a_i ) 362 ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 363 CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 364 & dta_bdy(jbdy)%h_i , dta_bdy(jbdy)%h_s , dta_bdy(jbdy)%a_i ) 365 ENDIF 359 CALL ice_var_itd( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 360 & dta_bdy(jbdy)%h_i , dta_bdy(jbdy)%h_s , dta_bdy(jbdy)%a_i ) 366 361 ENDIF 367 362 #endif -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/BDY/bdydyn2d.F90
r10529 r11422 187 187 ! Use characteristics method instead 188 188 zflag = ABS(flagu) 189 zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii m1,ij)189 zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij) 190 190 pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 191 191 END DO … … 205 205 ! Use characteristics method instead 206 206 zflag = ABS(flagv) 207 zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij m1)207 zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv)) 208 208 pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 209 209 END DO -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/BDY/bdyice.F90
r10909 r11422 57 57 INTEGER :: jbdy ! BDY set index 58 58 !!---------------------------------------------------------------------- 59 ! 60 IF( ln_timing ) CALL timing_start('bdy_ice_thd') 59 ! controls 60 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 61 IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 61 62 ! 62 63 CALL ice_var_glo2eqv … … 78 79 CALL ice_var_agg(1) 79 80 ! 80 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 81 IF( ln_timing ) CALL timing_stop('bdy_ice_thd') 81 ! controls 82 IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 83 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 84 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 82 85 ! 83 86 END SUBROUTINE bdy_ice … … 148 151 jpbound = 0 ; ib = ji ; jb = jj 149 152 ! 150 IF( u_ice(ji +1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 ; jb = jj151 IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji +1,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 ; jb = jj152 IF( v_ice(ji ,jj +1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; ib = ji; jb = jj+1153 IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj +1,1) == 0. ) jpbound = 1 ; ib = ji; jb = jj-1153 IF( u_ice(ji ,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 154 IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji ,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 155 IF( v_ice(ji ,jj ) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; jb = jj+1 156 IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj ,1) == 0. ) jpbound = 1 ; jb = jj-1 154 157 ! 155 158 IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DIA/diawri.F90
r10425 r11422 210 210 ENDIF 211 211 212 IF( ln_zad_Aimp ) wn = wn + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 213 ! 212 214 CALL iom_put( "woce", wn ) ! vertical velocity 213 215 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value … … 220 222 IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 221 223 ENDIF 224 ! 225 IF( ln_zad_Aimp ) wn = wn - wi ! Remove implicit part of vertical velocity that was added for diagnostic output 222 226 223 227 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. … … 842 846 CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress 843 847 844 CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current 848 IF( ln_zad_Aimp ) THEN 849 CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current 850 ELSE 851 CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current 852 ENDIF 845 853 CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef. 846 854 CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef. … … 903 911 CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity 904 912 CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity 905 CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity 913 IF( ln_zad_Aimp ) THEN 914 CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity 915 ELSE 916 CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity 917 ENDIF 906 918 IF( ALLOCATED(ahtu) ) THEN 907 919 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DOM/dommsk.F90
r10425 r11422 142 142 ENDIF 143 143 END DO 144 END DO 145 !SF add here lbc_lnk: bug not still understood : cause now domain configuration is read ! 146 !!gm I don't understand why... 144 END DO 145 ! 146 ! the following call is mandatory 147 ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...) 147 148 CALL lbc_lnk( 'dommsk', tmask , 'T', 1._wp ) ! Lateral boundary conditions 148 149 -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DOM/domvvl.F90
r11002 r11422 327 327 END DO 328 328 ! 329 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !330 ! ! ------baroclinic part------ !329 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 330 ! ! ------baroclinic part------ ! 331 331 ! I - initialization 332 332 ! ================== -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DYN/dynhpg.F90
r10491 r11422 37 37 USE trd_oce ! trends: ocean variables 38 38 USE trddyn ! trend manager: dynamics 39 !jcUSE zpshde ! partial step: hor. derivative (zps_hde routine)39 USE zpshde ! partial step: hor. derivative (zps_hde routine) 40 40 ! 41 41 USE in_out_manager ! I/O manager … … 338 338 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 339 339 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 340 REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv 340 341 !!---------------------------------------------------------------------- 341 342 ! … … 346 347 ENDIF 347 348 348 ! Partial steps: bottom beforehorizontal gradient of t, s, rd at the last ocean level349 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru ,grv )349 ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level 350 CALL zps_hde( kt, jpts, tsn, zgtsu, zgtsv, rhd, zgru , zgrv ) 350 351 351 352 ! Local constant initialization … … 385 386 END DO 386 387 387 ! partial steps correction at the last level (use gru &grv computed in zpshde.F90)388 ! partial steps correction at the last level (use zgru & zgrv computed in zpshde.F90) 388 389 DO jj = 2, jpjm1 389 390 DO ji = 2, jpim1 … … 395 396 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 396 397 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 397 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)398 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 398 399 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 399 400 ENDIF … … 401 402 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 402 403 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 403 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)404 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 404 405 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 405 406 ENDIF -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DYN/dynzdf.F90
r10364 r11422 170 170 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 171 171 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 172 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )173 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )172 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 173 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 174 174 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 175 175 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) … … 185 185 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 186 186 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 187 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )188 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )187 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 188 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 189 189 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 190 190 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) … … 199 199 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 200 200 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 201 zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) )201 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 202 202 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 203 203 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) … … 336 336 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 337 337 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 338 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )339 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)338 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 339 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 340 340 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 341 341 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) … … 351 351 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 352 352 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 353 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )354 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)353 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 354 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 355 355 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 356 356 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) … … 365 365 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 366 366 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 367 zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) )367 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 368 368 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 369 369 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/DYN/sshwzv.F90
r10907 r11422 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection 11 12 !!---------------------------------------------------------------------- 12 13 … … 279 280 !! : wi : now vertical velocity (for implicit treatment) 280 281 !! 281 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 282 !! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 283 !! implicit scheme for vertical advection in oceanic modeling. 284 !! Ocean Modelling, 91, 38-69. 282 285 !!---------------------------------------------------------------------- 283 286 INTEGER, INTENT(in) :: kt ! time step 284 287 ! 285 288 INTEGER :: ji, jj, jk ! dummy loop indices 286 REAL(wp) :: zCu, zcff, z1_e3 w! local scalars289 REAL(wp) :: zCu, zcff, z1_e3t ! local scalars 287 290 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 288 REAL(wp) , PARAMETER :: Cu_max = 0. 27! local parameters291 REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters 289 292 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 290 293 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters … … 297 300 IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 298 301 IF(lwp) WRITE(numout,*) '~~~~~ ' 299 ENDIF 300 ! 301 ! 302 DO jk = 1, jpkm1 ! calculate Courant numbers 303 DO jj = 2, jpjm1 304 DO ji = 2, fs_jpim1 ! vector opt. 305 z1_e3w = 1._wp / e3w_n(ji,jj,jk) 306 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & ! 2*rdt and not r2dt (for restartability) 307 & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 308 & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 309 & * r1_e1e2t(ji,jj) & 310 & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 311 & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 312 & * r1_e1e2t(ji,jj) & 313 & ) * z1_e3w 302 wi(:,:,:) = 0._wp 303 ENDIF 304 ! 305 ! Calculate Courant numbers 306 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 307 DO jk = 1, jpkm1 308 DO jj = 2, jpjm1 309 DO ji = 2, fs_jpim1 ! vector opt. 310 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 311 ! 2*rdt and not r2dt (for restartability) 312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 313 & + ( MAX( e2u(ji ,jj)*e3u_n(ji ,jj,jk)*un(ji ,jj,jk) + un_td(ji ,jj,jk), 0._wp ) - & 314 & MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) ) & 315 & * r1_e1e2t(ji,jj) & 316 & + ( MAX( e1v(ji,jj )*e3v_n(ji,jj ,jk)*vn(ji,jj ,jk) + vn_td(ji,jj ,jk), 0._wp ) - & 317 & MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) ) & 318 & * r1_e1e2t(ji,jj) & 319 & ) * z1_e3t 320 END DO 314 321 END DO 315 322 END DO 316 END DO 323 ELSE 324 DO jk = 1, jpkm1 325 DO jj = 2, jpjm1 326 DO ji = 2, fs_jpim1 ! vector opt. 327 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 328 ! 2*rdt and not r2dt (for restartability) 329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 330 & + ( MAX( e2u(ji ,jj)*e3u_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 331 & MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 332 & * r1_e1e2t(ji,jj) & 333 & + ( MAX( e1v(ji,jj )*e3v_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 334 & MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 335 & * r1_e1e2t(ji,jj) & 336 & ) * z1_e3t 337 END DO 338 END DO 339 END DO 340 ENDIF 317 341 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 318 342 ! … … 320 344 ! 321 345 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 322 DO jk = 1, jpkm1! or scan Courant criterion and partition346 DO jk = jpkm1, 2, -1 ! or scan Courant criterion and partition 323 347 DO jj = 1, jpj ! w where necessary 324 348 DO ji = 1, jpi 325 349 ! 326 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 350 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 351 ! alt: 352 ! IF ( wn(ji,jj,jk) > 0._wp ) THEN 353 ! zCu = Cu_adv(ji,jj,jk) 354 ! ELSE 355 ! zCu = Cu_adv(ji,jj,jk-1) 356 ! ENDIF 327 357 ! 328 358 IF( zCu <= Cu_min ) THEN !<-- Fully explicit … … 339 369 wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 340 370 ! 341 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient 371 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl 342 372 END DO 343 373 END DO 344 374 END DO 375 Cu_adv(:,:,1) = 0._wp 345 376 ELSE 346 377 ! Fully explicit everywhere 347 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient 378 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl 348 379 wi (:,:,:) = 0._wp 349 380 ENDIF -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/LBC/mppini.F90
r10615 r11422 669 669 ! 670 670 CALL mpp_init_ioipsl ! Prepare NetCDF output file (if necessary) 671 ! 672 IF ( ln_nnogather) THEN671 ! 672 IF (( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ).AND.( ln_nnogather )) THEN 673 673 CALL mpp_init_nfdcom ! northfold neighbour lists 674 674 IF (llwrtlay) THEN -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/SBC/sbcapr.F90
r10425 r11422 26 26 PUBLIC sbc_apr_init ! routine called in sbcmod 27 27 28 ! !!* namsbc_apr namelist (Atmospheric PRessure) *29 LOGICAL, PUBLIC :: ln_apr_obc !: inverse barometer added to OBC ssh data30 LOGICAL, PUBLIC :: ln_ref_apr !: ref. pressure: global mean Patm (F) or a constant (F)31 REAL(wp) :: rn_pref ! reference atmospheric pressure [N/m2]28 ! !!* namsbc_apr namelist (Atmospheric PRessure) * 29 LOGICAL, PUBLIC :: ln_apr_obc = .false. !: inverse barometer added to OBC ssh data 30 LOGICAL, PUBLIC :: ln_ref_apr !: ref. pressure: global mean Patm (F) or a constant (F) 31 REAL(wp) :: rn_pref ! reference atmospheric pressure [N/m2] 32 32 33 33 REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) :: ssh_ib ! Inverse barometer now sea surface height [m] -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/STO/stopar.F90
r10425 r11422 270 270 IF(lwm) WRITE ( numond, namsto ) 271 271 272 IF( .NOT.ln_ rststo) THEN ! no use of stochastic parameterization272 IF( .NOT.ln_sto_eos ) THEN ! no use of stochastic parameterization 273 273 IF(lwp) THEN 274 274 WRITE(numout,*) -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/TRA/traadv_fct.F90
r10425 r11422 21 21 USE diaar5 ! AR5 diagnostics 22 22 USE phycst , ONLY : rau0_rcp 23 USE zdf_oce , ONLY : ln_zad_Aimp 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 86 87 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 87 88 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 89 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup 90 LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection 88 91 !!---------------------------------------------------------------------- 89 92 ! … … 97 100 l_hst = .FALSE. 98 101 l_ptr = .FALSE. 102 ll_zAimp = .FALSE. 99 103 IF( ( cdtype =='TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 100 104 IF( cdtype =='TRA' .AND. ln_diaptr ) l_ptr = .TRUE. … … 116 120 ! 117 121 zwi(:,:,:) = 0._wp 122 ! 123 ! If adaptive vertical advection, check if it is needed on this PE at this time 124 IF( ln_zad_Aimp ) THEN 125 IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 126 END IF 127 ! If active adaptive vertical advection, build tridiagonal matrix 128 IF( ll_zAimp ) THEN 129 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 130 DO jk = 1, jpkm1 131 DO jj = 2, jpjm1 132 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 133 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk) 134 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t_a(ji,jj,jk) 135 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk) 136 END DO 137 END DO 138 END DO 139 END IF 118 140 ! 119 141 DO jn = 1, kjpt !== loop over the tracers ==! … … 169 191 END DO 170 192 END DO 193 194 IF ( ll_zAimp ) THEN 195 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 196 ! 197 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 198 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 199 DO jj = 2, jpjm1 200 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 202 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 203 ztw(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 204 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 205 END DO 206 END DO 207 END DO 208 DO jk = 1, jpkm1 209 DO jj = 2, jpjm1 210 DO ji = fs_2, fs_jpim1 ! vector opt. 211 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 212 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 213 END DO 214 END DO 215 END DO 216 ! 217 END IF 171 218 ! 172 219 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) … … 277 324 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 278 325 ENDIF 326 ! 327 IF ( ll_zAimp ) THEN 328 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 329 DO jj = 2, jpjm1 330 DO ji = fs_2, fs_jpim1 ! vector opt. 331 ! ! total intermediate advective trends 332 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 333 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 334 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 335 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 336 END DO 337 END DO 338 END DO 339 ! 340 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 341 ! 342 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 343 DO jj = 2, jpjm1 344 DO ji = fs_2, fs_jpim1 ! vector opt. 345 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 346 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 347 zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 348 END DO 349 END DO 350 END DO 351 END IF 279 352 ! 280 353 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. ) … … 289 362 DO jj = 2, jpjm1 290 363 DO ji = fs_2, fs_jpim1 ! vector opt. 291 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 292 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 293 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 294 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 295 END DO 296 END DO 297 END DO 364 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 365 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 366 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 367 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 368 zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 369 END DO 370 END DO 371 END DO 372 ! 373 IF ( ll_zAimp ) THEN 374 ! 375 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 376 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 377 DO jj = 2, jpjm1 378 DO ji = fs_2, fs_jpim1 ! vector opt. 379 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 380 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 381 ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 382 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 383 END DO 384 END DO 385 END DO 386 DO jk = 1, jpkm1 387 DO jj = 2, jpjm1 388 DO ji = fs_2, fs_jpim1 ! vector opt. 389 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 390 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 391 END DO 392 END DO 393 END DO 394 END IF 298 395 ! 299 396 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport … … 318 415 END DO ! end of tracer loop 319 416 ! 417 IF ( ll_zAimp ) THEN 418 DEALLOCATE( zwdia, zwinf, zwsup ) 419 ENDIF 320 420 IF( l_trd .OR. l_hst ) THEN 321 421 DEALLOCATE( ztrdx, ztrdy, ztrdz ) -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/TRA/traqsr.F90
r10425 r11422 168 168 DO jj = 2, jpjm1 ! Separation in R-G-B depending of the surface Chl 169 169 DO ji = fs_2, fs_jpim1 170 zchl = sf_chl(1)%fnow(ji,jj,1)170 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 171 171 zCtot = 40.6 * zchl**0.459 172 172 zze = 568.2 * zCtot**(-0.746) -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/step.F90
r10364 r11422 165 165 CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation 166 166 167 !!jc: fs simplification168 !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)169 !! but ensures reproductible results170 !! with previous versions using split-explicit free surface171 IF( ln_zps .AND. .NOT. ln_isfcav ) &172 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient173 & rhd, gru , grv ) ! of t, s, rd at the last ocean level174 IF( ln_zps .AND. ln_isfcav ) &175 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)176 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level177 !!jc: fs simplification178 167 179 168 ua(:,:,:) = 0._wp ! set dynamics trends to zero … … 198 187 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 199 188 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 189 ENDIF 190 CALL dyn_zdf ( kstp ) ! vertical diffusion 191 192 IF( ln_dynspg_ts ) THEN 200 193 CALL wzv ( kstp ) ! now cross-level velocity 201 194 IF( ln_zad_Aimp ) CALL wAimp ( kstp ) ! Adaptive-implicit vertical advection partitioning 202 195 ENDIF 203 204 CALL dyn_zdf ( kstp ) ! vertical diffusion205 196 206 197 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
NEMO/branches/2019/fix_vvl_ticket1791/src/OCE/stpctl.F90
r10570 r11422 96 96 IF( ln_zad_Aimp ) THEN 97 97 istatus = NF90_DEF_VAR( idrun, 'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 98 istatus = NF90_DEF_VAR( idrun, 'C u_max', NF90_DOUBLE, (/ idtime /), idc1 )98 istatus = NF90_DEF_VAR( idrun, 'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 ) 99 99 ENDIF 100 100 istatus = NF90_ENDDEF(idrun) … … 123 123 IF( ln_zad_Aimp ) THEN 124 124 zmax(8) = MAXVAL( ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 125 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! cell Courant no. max125 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 126 126 ENDIF 127 127 !
Note: See TracChangeset
for help on using the changeset viewer.
