Changeset 12489 for NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
- Timestamp:
- 2020-02-28T16:55:11+01:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
r12377 r12489 72 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 73 73 ! 74 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_ baro <= 2.5 nn_baro75 REAL(wp),SAVE :: r dtbt! Barotropic time step74 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e 75 REAL(wp),SAVE :: rDt_e ! Barotropic time step 76 76 ! 77 77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 102 102 ierr(:) = 0 103 103 ! 104 ALLOCATE( wgtbtp1(3*nn_ baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )104 ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) ) 105 105 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 106 106 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) … … 150 150 LOGICAL :: ll_init ! =T : special startup of 2d equations 151 151 INTEGER :: noffset ! local integers : time offset for bdy update 152 REAL(wp) :: r1_ 2dt_b, z1_hu, z1_hv ! local scalars152 REAL(wp) :: r1_Dt_b, z1_hu, z1_hv ! local scalars 153 153 REAL(wp) :: za0, za1, za2, za3 ! - - 154 154 REAL(wp) :: zztmp, zldg ! - - … … 180 180 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 181 181 ! ! inverse of baroclinic time step 182 IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt ) 183 ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt ) 184 ENDIF 182 r1_Dt_b = 1._wp / rDt 185 183 ! 186 184 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 185 ll_fw_start = .FALSE. 188 186 ! ! time offset in steps for bdy data update 189 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_ baro187 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 190 188 ELSE ; noffset = 0 191 189 ENDIF … … 198 196 IF(lwp) WRITE(numout,*) 199 197 ! 200 IF( neuler == 0) ll_init=.TRUE.201 ! 202 IF( ln_bt_fw .OR. neuler == 0) THEN198 IF( l_1st_euler ) ll_init=.TRUE. 199 ! 200 IF( ln_bt_fw .OR. l_1st_euler ) THEN 203 201 ll_fw_start =.TRUE. 204 202 noffset = 0 … … 209 207 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 210 208 ! 211 ENDIF 212 ! 213 ! If forward start at previous time step, and centered integration, 214 ! then update averaging weights: 215 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 216 ll_fw_start=.FALSE. 217 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 218 ENDIF 219 ! 220 209 ELSEIF( kt == nit000 + 1 ) THEN !* initialisation 2nd time-step 210 ! 211 IF( .NOT.ln_bt_fw ) THEN 212 ! If we did an Euler timestep on the first timestep we need to reset ll_fw_start 213 ! and the averaging weights. We don't have an easy way of telling whether we did 214 ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false. 215 ! at the end of the first timestep) so just do this in all cases. 216 ll_fw_start = .FALSE. 217 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 218 ENDIF 219 ! 220 ENDIF 221 ! 221 222 ! ----------------------------------------------------------------------------- 222 223 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 302 303 IF( ln_bt_fw ) THEN ! Add wind forcing 303 304 DO_2D_00_00 304 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_r au0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)305 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_r au0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)305 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 306 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) 306 307 END_2D 307 308 ELSE 308 zztmp = r1_r au0 * r1_2309 zztmp = r1_rho0 * r1_2 309 310 DO_2D_00_00 310 311 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) … … 319 320 ! ! --------------------------------------------------- ! 320 321 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 321 zssh_frc(:,:) = r1_r au0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )322 zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 322 323 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 323 zztmp = r1_r au0 * r1_2324 zztmp = r1_rho0 * r1_2 324 325 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 325 326 & - rnf(:,:) - rnf_b(:,:) & … … 428 429 ! Update tide potential at the beginning of current time substep 429 430 IF( ln_tide_pot .AND. ln_tide ) THEN 430 zt0substep = REAL(nsec_day, wp) - 0.5_wp*r dt + (jn + noffset - 1) * rdt / REAL(nn_baro, wp)431 zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_e, wp) 431 432 CALL upd_tide(zt0substep, Kmm) 432 433 END IF … … 494 495 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV ) 495 496 #endif 496 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, r dtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV497 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 497 498 498 499 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where … … 509 510 DO_2D_00_00 510 511 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 511 ssha_e(ji,jj) = ( sshn_e(ji,jj) - r dtbt* ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)512 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 512 513 END_2D 513 514 ! … … 599 600 DO_2D_00_00 600 601 ua_e(ji,jj) = ( un_e(ji,jj) & 601 & + r dtbt* ( zu_spg(ji,jj) &602 & + rDt_e * ( zu_spg(ji,jj) & 602 603 & + zu_trd(ji,jj) & 603 604 & + zu_frc(ji,jj) ) & … … 605 606 606 607 va_e(ji,jj) = ( vn_e(ji,jj) & 607 & + r dtbt* ( zv_spg(ji,jj) &608 & + rDt_e * ( zv_spg(ji,jj) & 608 609 & + zv_trd(ji,jj) & 609 610 & + zv_frc(ji,jj) ) & … … 624 625 ! 625 626 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 626 & + r dtbt* ( zhu_bck * zu_spg (ji,jj) & !627 & + rDt_e * ( zhu_bck * zu_spg (ji,jj) & ! 627 628 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 628 629 & + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu 629 630 ! 630 631 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 631 & + r dtbt* ( zhv_bck * zv_spg (ji,jj) & !632 & + rDt_e * ( zhv_bck * zv_spg (ji,jj) & ! 632 633 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 633 634 & + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv … … 637 638 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 638 639 DO_2D_00_00 639 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - r dtbt* zCdU_u(ji,jj) * hur_e(ji,jj))640 va_e(ji,jj) = va_e(ji,jj) /(1.0 - r dtbt* zCdU_v(ji,jj) * hvr_e(ji,jj))640 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 641 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 641 642 END_2D 642 643 ENDIF … … 701 702 ! Set advection velocity correction: 702 703 IF (ln_bt_fw) THEN 703 IF( .NOT.( kt == nit000 .AND. neuler==0) ) THEN704 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 704 705 DO_2D_11_11 705 706 zun_save = un_adv(ji,jj) 706 707 zvn_save = vn_adv(ji,jj) 707 708 ! ! apply the previously computed correction 708 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )709 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )709 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) ) 710 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) ) 710 711 ! ! Update corrective fluxes for next time step 711 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )712 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )712 un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 713 vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 713 714 ! ! Save integrated transport for next computation 714 715 ub2_b(ji,jj) = zun_save … … 728 729 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 729 730 DO jk=1,jpkm1 730 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_ 2dt_b731 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_ 2dt_b731 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b 732 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b 732 733 END DO 733 734 ELSE … … 744 745 ! 745 746 DO jk=1,jpkm1 746 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_ 2dt_b747 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_ 2dt_b747 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 748 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 748 749 END DO 749 750 ! Save barotropic velocities not transport: … … 808 809 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 809 810 INTEGER, INTENT(inout) :: jpit ! cycle length 810 REAL(wp), DIMENSION(3*nn_ baro), INTENT(inout) :: zwgt1, & ! Primary weights811 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 811 812 zwgt2 ! Secondary weights 812 813 … … 820 821 ! Set time index when averaged value is requested 821 822 IF (ll_fw) THEN 822 jic = nn_ baro823 jic = nn_e 823 824 ELSE 824 jic = 2 * nn_ baro825 jic = 2 * nn_e 825 826 ENDIF 826 827 … … 828 829 IF (ll_av) THEN 829 830 ! Define simple boxcar window for primary weights 830 ! (width = nn_ baro, centered around jic)831 ! (width = nn_e, centered around jic) 831 832 SELECT CASE ( nn_bt_flt ) 832 833 CASE( 0 ) ! No averaging … … 834 835 jpit = jic 835 836 836 CASE( 1 ) ! Boxcar, width = nn_ baro837 DO jn = 1, 3*nn_ baro838 za1 = ABS(float(jn-jic))/float(nn_ baro)837 CASE( 1 ) ! Boxcar, width = nn_e 838 DO jn = 1, 3*nn_e 839 za1 = ABS(float(jn-jic))/float(nn_e) 839 840 IF (za1 < 0.5_wp) THEN 840 841 zwgt1(jn) = 1._wp … … 843 844 ENDDO 844 845 845 CASE( 2 ) ! Boxcar, width = 2 * nn_ baro846 DO jn = 1, 3*nn_ baro847 za1 = ABS(float(jn-jic))/float(nn_ baro)846 CASE( 2 ) ! Boxcar, width = 2 * nn_e 847 DO jn = 1, 3*nn_e 848 za1 = ABS(float(jn-jic))/float(nn_e) 848 849 IF (za1 < 1._wp) THEN 849 850 zwgt1(jn) = 1._wp … … 889 890 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 890 891 ! ! --------------- 891 IF( ln_rstart .AND. ln_bt_fw .AND. ( neuler/=0) ) THEN !* Read the restart file892 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 892 893 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios ) 893 894 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios ) … … 975 976 976 977 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 977 IF( ln_bt_auto ) nn_ baro = CEILING( rdt / rn_bt_cmax * zcmax)978 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 978 979 979 r dtbt = rdt / REAL( nn_baro, wp )980 zcmax = zcmax * r dtbt980 rDt_e = rn_Dt / REAL( nn_e , wp ) 981 zcmax = zcmax * rDt_e 981 982 ! Print results 982 983 IF(lwp) WRITE(numout,*) … … 984 985 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 985 986 IF( ln_bt_auto ) THEN 986 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_ baro'987 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e ' 987 988 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 988 989 ELSE 989 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_ baro in namelist nn_baro = ', nn_baro990 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 990 991 ENDIF 991 992 992 993 IF(ln_bt_av) THEN 993 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_ barotime steps is on '994 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' 994 995 ELSE 995 996 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' … … 1011 1012 SELECT CASE ( nn_bt_flt ) 1012 1013 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1013 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_ baro'1014 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_ baro'1014 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1015 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1015 1016 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1016 1017 END SELECT 1017 1018 ! 1018 1019 IF(lwp) WRITE(numout,*) ' ' 1019 IF(lwp) WRITE(numout,*) ' nn_ baro = ', nn_baro1020 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', r dtbt1020 IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e 1021 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e 1021 1022 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1022 1023 ! … … 1030 1031 ENDIF 1031 1032 IF( zcmax>0.9_wp ) THEN 1032 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_ baro!' )1033 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1033 1034 ENDIF 1034 1035 ! … … 1429 1430 ! 1430 1431 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1431 zztmp = -1._wp / r dtbt1432 zztmp = -1._wp / rDt_e 1432 1433 DO_2D_00_00 1433 1434 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
Note: See TracChangeset
for help on using the changeset viewer.
