Changeset 14063 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
- Timestamp:
- 2020-12-03T17:45:05+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
r14017 r14063 306 306 ENDIF 307 307 ! 308 ! != Add atmospheric pressureforcing =!309 ! ! ------------------ ----------------!310 IF( ln_bt_fw ) THEN ! Add wind forcing308 ! != Add wind forcing =! 309 ! ! ------------------ ! 310 IF( ln_bt_fw ) THEN 311 311 DO_2D( 0, 0, 0, 0 ) 312 312 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) … … 386 386 ! 387 387 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 388 zhup2_e(:,:) = hu (:,:,Kmm)389 zhvp2_e(:,:) = hv (:,:,Kmm)390 zhtp2_e(:,:) = ht (:,:)391 ENDIF 392 ! 393 IF (ln_bt_fw) THEN! FORWARD integration: start from NOW fields394 sshn_e(:,:) = pssh (:,:,Kmm)388 zhup2_e(:,:) = hu_0(:,:) 389 zhvp2_e(:,:) = hv_0(:,:) 390 zhtp2_e(:,:) = ht_0(:,:) 391 ENDIF 392 ! 393 IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields 394 sshn_e(:,:) = pssh (:,:,Kmm) 395 395 un_e (:,:) = puu_b(:,:,Kmm) 396 396 vn_e (:,:) = pvv_b(:,:,Kmm) … … 401 401 hvr_e (:,:) = r1_hv(:,:,Kmm) 402 402 ELSE ! CENTRED integration: start from BEFORE fields 403 sshn_e(:,:) = pssh (:,:,Kbb)403 sshn_e(:,:) = pssh (:,:,Kbb) 404 404 un_e (:,:) = puu_b(:,:,Kbb) 405 405 vn_e (:,:) = pvv_b(:,:,Kbb) … … 412 412 ! 413 413 ! Initialize sums: 414 puu_b 415 pvv_b 414 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 415 pvv_b (:,:,Kaa) = 0._wp 416 416 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level 417 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop418 vn_adv(:,:) = 0._wp417 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 418 vn_adv(:,:) = 0._wp 419 419 ! 420 420 IF( ln_wd_dl ) THEN … … 475 475 ! 476 476 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 477 #if defined key_qcoTest_FluxForm 478 ! ! 'key_qcoTest_FluxForm' : simple ssh average 479 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 480 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj) 481 END_2D 482 DO_2D( 1, 0, 1, 1 ) 483 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj) 484 END_2D 485 #else 486 ! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 477 487 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 478 488 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & … … 485 495 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 486 496 END_2D 497 #endif 487 498 ! 488 499 ENDIF … … 541 552 ! Sea Surface Height at u-,v-points (vvl case only) 542 553 IF( .NOT.ln_linssh ) THEN 554 #if defined key_qcoTest_FluxForm 555 ! ! 'key_qcoTest_FluxForm' : simple ssh average 556 DO_2D( 1, 1, 1, 0 ) 557 zsshu_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji+1,jj ) ) * ssumask(ji,jj) 558 END_2D 559 DO_2D( 1, 0, 1, 1 ) 560 zsshv_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji ,jj+1) ) * ssvmask(ji,jj) 561 END_2D 562 #else 543 563 DO_2D( 0, 0, 0, 0 ) 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 546 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 547 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 548 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 549 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 550 END_2D 564 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 565 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) * ssumask(ji,jj) 566 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 567 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) * ssvmask(ji,jj) 568 END_2D 569 #endif 551 570 ENDIF 552 571 ! … … 624 643 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 625 644 ! ! backward interpolated depth used in spg terms at jn+1/2 645 #if defined key_qcoTest_FluxForm 646 ! ! 'key_qcoTest_FluxForm' : simple ssh average 647 zhu_bck = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj) 648 zhv_bck = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj) 649 #else 626 650 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 627 651 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 628 652 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 629 653 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 654 #endif 630 655 ! ! inverse depth at jn+1 631 656 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) … … 646 671 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 647 672 DO_2D( 0, 0, 0, 0 ) 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))673 ua_e(ji,jj) = ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) ) 674 va_e(ji,jj) = va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) ) 650 675 END_2D 651 676 ENDIF 652 677 653 678 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)655 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1))656 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)657 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1))679 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 680 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 681 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 682 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 658 683 ENDIF 659 684 ! … … 743 768 ELSE 744 769 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 770 #if defined key_qcoTest_FluxForm 771 ! ! 'key_qcoTest_FluxForm' : simple ssh average 745 772 DO_2D( 1, 0, 1, 0 ) 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 748 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 749 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 750 & * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 751 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 752 END_2D 773 zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj) 774 zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj) 775 END_2D 776 #else 777 DO_2D( 1, 0, 1, 0 ) 778 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 779 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 780 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 781 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 782 END_2D 783 #endif 753 784 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 754 785 ! 755 786 DO jk=1,jpkm1 756 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 757 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 787 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) & 788 & * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 789 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) & 790 & * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 758 791 END DO 759 792 ! Save barotropic velocities not transport: … … 899 932 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 900 933 ! ! --------------- 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler)) THEN !* Read the restart file934 IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file 902 935 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp ) 903 936 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp ) … … 1060 1093 !! although they should be updated in the variable volume case. Not a big approximation. 1061 1094 !! To remove this approximation, copy lines below inside barotropic loop 1062 !! and update depths at T- F points (ht and zhf resp.) at each barotropic time step1095 !! and update depths at T- points (ht) at each barotropic time step 1063 1096 !! 1064 1097 !! Compute zwz = f / ( height of the water colomn ) … … 1067 1100 INTEGER :: ji ,jj, jk ! dummy loop indices 1068 1101 REAL(wp) :: z1_ht 1069 REAL(wp), DIMENSION(jpi,jpj) :: zhf1070 1102 !!---------------------------------------------------------------------- 1071 1103 ! 1072 1104 SELECT CASE( nvor_scheme ) 1073 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme1074 SELECT CASE( nn_e en_e3f )!* ff_f/e3 at F-point1105 CASE( np_EEN, np_ENE, np_ENS , np_MIX ) != schemes using the same e3f definition 1106 SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point 1075 1107 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1076 DO_2D( 1, 0, 1, 0 )1077 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) +&1078 & ht(ji ,jj ) + ht(ji+1,jj )) * 0.25_wp1108 DO_2D( 0, 0, 0, 0 ) 1109 zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) & 1110 & + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1079 1111 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1080 1112 END_2D 1081 1113 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1082 DO_2D( 1, 0, 1, 0 )1083 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) &1084 & + ht (ji ,jj ) + ht(ji+1,jj ) ) &1085 & / ( MAX( 1._wp, ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &1086 & + ssmask(ji ,jj ) + ssmask(ji+1,jj )) )1114 DO_2D( 0, 0, 0, 0 ) 1115 zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) & 1116 & + ht(ji,jj ) + ht(ji+1,jj ) ) & 1117 & / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) & 1118 & + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) ) 1087 1119 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1088 1120 END_2D 1089 1121 END SELECT 1090 1122 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1091 ! 1092 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1123 END SELECT 1124 ! 1125 SELECT CASE( nvor_scheme ) 1126 CASE( np_EEN ) 1127 ! 1128 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1093 1129 DO_2D( 0, 1, 0, 1 ) 1094 1130 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) … … 1098 1134 END_2D 1099 1135 ! 1100 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme1101 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ;ftsw(1,:) = 0._wp1136 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1137 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1102 1138 DO_2D( 0, 1, 0, 1 ) 1103 1139 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) … … 1108 1144 END_2D 1109 1145 ! 1110 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !1111 !1112 zwz(:,:) = 0._wp1113 zhf(:,:) = 0._wp1114 1115 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed1116 !!gm A priori a better value should be something like :1117 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)1118 !!gm divided by the sum of the corresponding mask1119 !!gm1120 !!1121 IF( .NOT.ln_sco ) THEN1122 1123 !!gm agree the JC comment : this should be done in a much clear way1124 1125 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case1126 ! Set it to zero for the time being1127 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level1128 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth1129 ! ENDIF1130 ! zhf(:,:) = gdepw_0(:,:,jk+1)1131 !1132 ELSE1133 !1134 !zhf(:,:) = hbatf(:,:)1135 DO_2D( 1, 0, 1, 0 )1136 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &1137 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &1138 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &1139 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )1140 END_2D1141 ENDIF1142 !1143 DO jj = 1, jpjm11144 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))1145 END DO1146 !1147 DO jk = 1, jpkm11148 DO jj = 1, jpjm11149 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)1150 END DO1151 END DO1152 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )1153 ! JC: TBC. hf should be greater than 01154 DO_2D( 1, 1, 1, 1 )1155 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)1156 END_2D1157 zwz(:,:) = ff_f(:,:) * zwz(:,:)1158 1146 END SELECT 1159 1147 1160 1148 END SUBROUTINE dyn_cor_2d_init 1161 1162 1149 1163 1150 … … 1352 1339 1353 1340 END SUBROUTINE wad_spg 1354 1355 1341 1356 1342
Note: See TracChangeset
for help on using the changeset viewer.
