- Timestamp:
- 2018-06-30T12:51:02+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90
r9598 r9863 56 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 57 57 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 58 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors59 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors58 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n ! baroclinic scale factors 59 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a ! baroclinic scale factors 60 60 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 61 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence … … 76 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 ALLOCATE( t ilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &79 & dt ilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &78 ALLOCATE( te3t_b(jpi,jpj,jpk) , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) , & 79 & dte3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , & 80 80 & STAT = dom_vvl_alloc ) 81 81 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) … … 103 103 !! - interpolate scale factors 104 104 !! 105 !! ** Action : - e3t_(n/b) and t ilde_e3t_(n/b)105 !! ** Action : - e3t_(n/b) and te3t_(n/b) 106 106 !! - Regrid: e3(u/v)_n 107 107 !! e3(u/v)_b … … 129 129 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 130 130 ! 131 ! ! Read or initialize e3t_(b/n), t ilde_e3t_(b/n) and hdiv_lf131 ! ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf 132 132 CALL dom_vvl_rst( nit000, 'READ' ) 133 133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all … … 280 280 !! 281 281 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 282 !! - t ilde_e3t_a: after increment of vertical scale factor282 !! - te3t_a: after increment of vertical scale factor 283 283 !! in z_tilde case 284 284 !! - e3(t/u/v)_a … … 353 353 ! II - after z_tilde increments of vertical scale factors 354 354 ! ======================================================= 355 t ilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms355 te3t_a(:,:,:) = 0._wp ! te3t_a used to store tendency terms 356 356 357 357 ! 1 - High frequency divergence term … … 359 359 IF( ln_vvl_ztilde ) THEN ! z_tilde case 360 360 DO jk = 1, jpkm1 361 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )361 te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 362 362 END DO 363 363 ELSE ! layer case 364 364 DO jk = 1, jpkm1 365 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)365 te3t_a(:,:,jk) = te3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 366 366 END DO 367 367 ENDIF … … 371 371 IF( ln_vvl_ztilde ) THEN 372 372 DO jk = 1, jpk 373 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)373 te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 374 374 END DO 375 375 ENDIF … … 383 383 DO ji = 1, fs_jpim1 ! vector opt. 384 384 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 385 & * ( t ilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )385 & * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj ,jk) ) 386 386 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 387 & * ( t ilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )387 & * ( te3t_b(ji,jj,jk) - te3t_b(ji ,jj+1,jk) ) 388 388 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 389 389 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 400 400 DO jj = 2, jpjm1 401 401 DO ji = fs_2, fs_jpim1 ! vector opt. 402 t ilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &402 te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 403 403 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 404 404 & ) * r1_e1e2t(ji,jj) … … 414 414 ! Leapfrog time stepping 415 415 ! ~~~~~~~~~~~~~~~~~~~~~~ 416 IF( neuler == 0 .AND. kt == nit000 ) THEN 417 z2dt = rdt 418 ELSE 419 z2dt = 2.0_wp * rdt 420 ENDIF 421 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 422 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 416 CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 417 te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 423 418 424 419 ! Maximum deformation control … … 426 421 ze3t(:,:,jpk) = 0._wp 427 422 DO jk = 1, jpkm1 428 ze3t(:,:,jk) = t ilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)423 ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 429 424 END DO 430 425 z_tmax = MAXVAL( ze3t(:,:,:) ) … … 446 441 ENDIF 447 442 IF (lwp) THEN 448 WRITE(numout, *) 'MAX( t ilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax443 WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 449 444 WRITE(numout, *) 'at i, j, k=', ijk_max 450 WRITE(numout, *) 'MIN( t ilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin445 WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 451 446 WRITE(numout, *) 'at i, j, k=', ijk_min 452 CALL ctl_warn('MAX( ABS( t ilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')447 CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 453 448 ENDIF 454 449 ENDIF 455 450 ! - ML - end test 456 451 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 457 t ilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )458 t ilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )452 te3t_a(:,:,:) = MIN( te3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 453 te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 459 454 460 455 ! … … 462 457 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 463 458 DO jk = 1, jpkm1 464 dt ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)459 dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 465 460 END DO 466 461 ! III - Barotropic repartition of the sea surface height over the baroclinic profile … … 470 465 ! i.e. locally and not spread over the water column. 471 466 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 472 zht(:,:) = 0. 467 zht(:,:) = 0._wp 473 468 DO jk = 1, jpkm1 474 zht(:,:) = zht(:,:) + t ilde_e3t_a(:,:,jk) * tmask(:,:,jk)469 zht(:,:) = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 475 470 END DO 476 471 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 477 472 DO jk = 1, jpkm1 478 dt ilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)473 dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 479 474 END DO 480 475 … … 484 479 ! ! ---baroclinic part--------- ! 485 480 DO jk = 1, jpkm1 486 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dt ilde_e3t_a(:,:,jk) * tmask(:,:,jk)481 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 487 482 END DO 488 483 ENDIF … … 494 489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 495 490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 496 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(t ilde_e3t_a))) =', z_tmax491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 497 492 END IF 498 493 ! … … 573 568 !! - recompute depths and water height fields 574 569 !! 575 !! ** Action : - e3t_(b/n), t ilde_e3t_(b/n) and e3(u/v)_n ready for next time step570 !! ** Action : - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step 576 571 !! - Recompute: 577 572 !! e3(u/v)_b … … 587 582 INTEGER, INTENT( in ) :: kt ! time step 588 583 ! 589 INTEGER :: ji, jj, jk ! dummy loop indices590 REAL(wp) :: zcoef 584 INTEGER :: ji, jj, jk ! dummy loop indices 585 REAL(wp) :: zcoef, ze3f ! local scalar 591 586 !!---------------------------------------------------------------------- 592 587 ! … … 605 600 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 606 601 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 607 IF( neuler == 0 .AND. kt == nit000) THEN608 t ilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)602 IF( l_1st_euler ) THEN 603 te3t_n(:,:,:) = te3t_a(:,:,:) 609 604 ELSE 610 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 611 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 605 DO jk = 1, jpk 606 DO jj = 1, jpj 607 DO ji = 1, jpi 608 ze3f = te3t_n(ji,jj,jk) & 609 & + atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 610 te3t_b(ji,jj,jk) = ze3f 611 te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 612 END DO 613 END DO 614 END DO 612 615 ENDIF 613 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)614 616 ENDIF 615 617 gdept_b(:,:,:) = gdept_n(:,:,:) … … 806 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios ) 807 809 ! 808 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. )809 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. )810 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. ) 811 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. ) 810 812 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 811 813 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 812 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. )814 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. ) 813 815 ! ! --------- ! 814 816 ! ! all cases ! … … 823 825 e3t_b(:,:,:) = e3t_0(:,:,:) 824 826 END WHERE 825 IF( neuler == 0) THEN827 IF( l_1st_euler ) THEN 826 828 e3t_b(:,:,:) = e3t_n(:,:,:) 827 829 ENDIF … … 829 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 830 832 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 831 IF(lwp) write(numout,*) ' neuler is forced to 0'833 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 832 834 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 833 835 e3t_n(:,:,:) = e3t_b(:,:,:) 834 neuler = 0836 l_1st_euler = .TRUE. 835 837 ELSE IF( id2 > 0 ) THEN 836 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 837 839 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 838 IF(lwp) write(numout,*) ' neuler is forced to 0'840 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 839 841 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 840 842 e3t_b(:,:,:) = e3t_n(:,:,:) 841 neuler = 0843 l_1st_euler = .TRUE. 842 844 ELSE 843 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 844 846 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 845 IF(lwp) write(numout,*) ' neuler is forced to 0'847 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 846 848 DO jk = 1, jpk 847 849 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & … … 850 852 END DO 851 853 e3t_b(:,:,:) = e3t_n(:,:,:) 852 neuler = 0854 l_1st_euler = .TRUE. 853 855 ENDIF 854 856 ! ! ----------- ! … … 862 864 ! ! ----------------------- ! 863 865 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 864 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', t ilde_e3t_b(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', t ilde_e3t_n(:,:,:), ldxios = lrxios )866 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios ) 867 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios ) 866 868 ELSE ! one at least array is missing 867 t ilde_e3t_b(:,:,:) = 0.0_wp868 t ilde_e3t_n(:,:,:) = 0.0_wp869 te3t_b(:,:,:) = 0.0_wp 870 te3t_n(:,:,:) = 0.0_wp 869 871 ENDIF 870 872 ! ! ------------ ! … … 942 944 943 945 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 944 t ilde_e3t_b(:,:,:) = 0._wp945 t ilde_e3t_n(:,:,:) = 0._wp946 te3t_b(:,:,:) = 0._wp 947 te3t_n(:,:,:) = 0._wp 946 948 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 947 949 END IF … … 960 962 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 961 963 ! ! ----------------------- ! 962 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', t ilde_e3t_b(:,:,:), ldxios = lwxios)963 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', t ilde_e3t_n(:,:,:), ldxios = lwxios)964 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 965 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 964 966 END IF 965 967 ! ! -------------!
Note: See TracChangeset
for help on using the changeset viewer.
