Changeset 5930 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2015-11-26T17:07:10+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5913 r5930 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 14 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts15 15 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts'split explicit free surface16 !! split explicit free surface 17 17 !!---------------------------------------------------------------------- 18 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables26 25 USE phycst ! physical constants 27 26 USE dynvor ! vorticity term 28 27 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 28 USE bdytides ! open boundary condition data 30 29 USE bdydyn2d ! open boundary conditions on barotropic variables 31 30 USE sbctide ! tides … … 68 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 68 70 ! Arrays below are saved to allow testing of the "no time averaging" option71 ! If this option is not retained, these could be replaced by temporary arrays72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays73 ubb_e, ub_e, &74 vbb_e, vb_e75 76 69 !! * Substitutions 77 70 # include "domzgr_substitute.h90" … … 88 81 !! *** routine dyn_spg_ts_alloc *** 89 82 !!---------------------------------------------------------------------- 90 INTEGER :: ierr( 3)83 INTEGER :: ierr(4) 91 84 !!---------------------------------------------------------------------- 92 85 ierr(:) = 0 93 86 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 97 91 98 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) … … 101 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 103 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 104 105 105 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc … … 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 150 REAL(wp) :: zu_spg, zv_spg ! - -151 REAL(wp) :: zhura, zhvra 152 REAL(wp) :: za0, za1, za2, za3 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 REAL(wp) :: zu_spg, zv_spg ! - - 151 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 155 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 156 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 158 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 164 164 ! !* Allocate temporary arrays 165 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)167 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 485 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 486 sshn_e(:,:) = sshn (:,:) 487 zun_e(:,:) = un_b (:,:)488 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 489 489 ! 490 490 hu_e (:,:) = hu (:,:) … … 494 494 ELSE ! CENTRED integration: start from BEFORE fields 495 495 sshn_e(:,:) = sshb (:,:) 496 zun_e(:,:) = ub_b (:,:)497 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 498 498 ! 499 499 hu_e (:,:) = hu_b (:,:) … … 509 509 va_b (:,:) = 0._wp 510 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 511 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop512 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 513 513 ! ! ==================== ! 514 514 DO jn = 1, icycle ! sub-time-step loop ! … … 534 534 535 535 ! Extrapolate barotropic velocities at step jit+0.5: 536 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)537 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)536 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 537 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 538 539 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 597 597 ! Sum over sub-time-steps to compute advective velocities 598 598 za2 = wgtbtp2(jn) 599 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)600 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 601 601 ! 602 602 ! Set next sea level: … … 733 733 ! 734 734 ! Add bottom stresses: 735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 737 ! 738 738 ! Surface pressure trend: … … 751 751 DO jj = 2, jpjm1 752 752 DO ji = fs_2, fs_jpim1 ! vector opt. 753 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 754 754 & + rdtbt * ( zwx(ji,jj) & 755 755 & + zu_trd(ji,jj) & … … 757 757 & ) * umask(ji,jj,1) 758 758 759 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 760 760 & + rdtbt * ( zwy(ji,jj) & 761 761 & + zv_trd(ji,jj) & … … 772 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 773 773 774 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 778 778 & ) * zhura 779 779 780 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 781 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 802 802 #if defined key_bdy 803 803 ! open boundaries 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 805 #endif 806 806 #if defined key_agrif … … 810 810 ! ! ---- 811 811 ubb_e (:,:) = ub_e (:,:) 812 ub_e (:,:) = zun_e(:,:)813 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 814 814 ! 815 815 vbb_e (:,:) = vb_e (:,:) 816 vb_e (:,:) = zvn_e(:,:)817 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 818 818 ! 819 819 sshbb_e(:,:) = sshb_e(:,:) … … 840 840 ! ----------------------------------------------------------------------------- 841 841 ! 842 ! At this stage ssha holds a time averaged value 843 ! ! Sea Surface Height at u-,v- and f-points 844 IF( lk_vvl ) THEN ! (required only in key_vvl case) 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 848 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 851 END IF 852 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 856 ENDIF 857 ! 858 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 860 DO jk=1,jpkm1 861 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 862 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 863 END DO 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 845 866 DO jj = 1, jpjm1 846 867 DO ji = 1, jpim1 ! NO Vector Opt. … … 854 875 END DO 855 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 856 ENDIF 857 ! 858 ! Set advection velocity correction: 859 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 860 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 861 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 862 ELSE 863 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 864 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 865 END IF 866 867 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 868 ub2_b(:,:) = zu_sum(:,:) 869 vb2_b(:,:) = zv_sum(:,:) 870 ENDIF 871 ! 872 ! Update barotropic trend: 873 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 874 DO jk=1,jpkm1 875 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 876 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 877 END DO 878 ELSE 877 ! 879 878 DO jk=1,jpkm1 880 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 915 914 ! 916 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 917 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)918 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 919 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 920 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1064 1063 ! 1065 1064 INTEGER :: ji ,jj 1066 INTEGER :: ios ! Local integer output status for namelist read1067 1065 REAL(wp) :: zxr2, zyr2, zcmax 1068 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1069 1067 !! 1070 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1071 & nn_baro, rn_bt_cmax, nn_bt_flt1072 1068 !!---------------------------------------------------------------------- 1073 1069 ! 1074 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1075 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1076 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1077 1078 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1079 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1080 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1081 IF(lwm) WRITE ( numond, namsplit ) 1082 ! 1083 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1084 1071 ! 1085 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1086 1073 ! 1087 IF (lk_vvl) THEN 1088 DO jj = 1, jpj 1089 DO ji =1, jpi 1090 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1091 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1092 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1093 END DO 1094 END DO 1095 ELSE 1096 DO jj = 1, jpj 1097 DO ji =1, jpi 1098 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 1101 END DO 1102 END DO 1103 ENDIF 1104 1074 DO jj = 1, jpj 1075 DO ji =1, jpi 1076 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1077 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1078 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1079 END DO 1080 END DO 1081 ! 1105 1082 zcmax = MAXVAL( zcu(:,:) ) 1106 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1107 1084 1108 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1109 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1086 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1110 1087 1111 1088 rdtbt = rdt / REAL( nn_baro , wp ) … … 1115 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1116 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1117 IF( ln_bt_ nn_auto ) THEN1118 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1094 IF( ln_bt_auto ) THEN 1095 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1119 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1120 1097 ELSE 1121 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1098 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1122 1099 ENDIF 1123 1100 … … 1164 1141 END SUBROUTINE dyn_spg_ts_init 1165 1142 1166 #else1167 !!---------------------------------------------------------------------------1168 !! Default case : Empty module No split explicit free surface1169 !!---------------------------------------------------------------------------1170 CONTAINS1171 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1172 dyn_spg_ts_alloc = 01173 END FUNCTION dyn_spg_ts_alloc1174 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1175 INTEGER, INTENT(in) :: kt1176 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1177 END SUBROUTINE dyn_spg_ts1178 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1179 INTEGER , INTENT(in) :: kt ! ocean time-step1180 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1181 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1182 END SUBROUTINE ts_rst1183 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1184 INTEGER , INTENT(in) :: kt ! ocean time-step1185 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1186 END SUBROUTINE dyn_spg_ts_init1187 #endif1188 1189 1143 !!====================================================================== 1190 1144 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
