- Timestamp:
- 2017-12-12T16:38:41+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8992 1 1 MODULE dynspg_ts 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out ! 4 2 5 !!====================================================================== 3 6 !! *** MODULE dynspg_ts *** … … 137 140 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 141 ! 139 LOGICAL :: ll_fw_start 140 LOGICAL :: ll_init 142 LOGICAL :: ll_fw_start ! if true, forward integration 143 LOGICAL :: ll_init ! if true, special startup of 2d equations 141 144 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 142 INTEGER :: ji, jj, jk, jn 143 INTEGER :: ikbu, ikbv, noffset 144 INTEGER :: iktu, iktv 145 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 INTEGER :: iktu, iktv ! local integers 145 148 REAL(wp) :: zmdi 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf 149 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 150 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 151 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 152 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 153 REAL(wp) :: zhura, zhvra ! - - 154 REAL(wp) :: za0, za1, za2, za3 ! - - 155 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 156 157 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 158 152 159 ! 153 160 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e … … 158 165 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 166 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 167 REAL(wp), POINTER, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 168 REAL(wp), POINTER, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 160 169 !!---------------------------------------------------------------------- 161 170 ! … … 169 178 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 170 179 CALL wrk_alloc( jpi,jpj, zhf ) 171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 180 IF( ln_wd_il ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 181 IF( ln_wd_dl ) CALL wrk_alloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2) 182 183 172 184 ! 173 185 zmdi=1.e+20 ! missing data indicator for masking … … 178 190 z1_2 = 0.5_wp 179 191 zraur = 1._wp / rau0 180 ! ! reciprocal of baroclinic time step 192 zwdramp = r_rn_wdmin1 ! simplest ramp 193 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 194 ! ! reciprocal of baroclinic time step 181 195 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 182 196 ELSE ; z2dt_bf = 2.0_wp * rdt … … 184 198 z1_2dt_b = 1.0_wp / z2dt_bf 185 199 ! 186 ll_init = ln_bt_av 200 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 201 ll_fw_start = .FALSE. 188 ! 202 ! ! time offset in steps for bdy data update 189 203 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 190 204 ELSE ; noffset = 0 191 205 ENDIF 192 206 ! 193 IF( kt == nit000 ) THEN 207 IF( kt == nit000 ) THEN !* initialisation 194 208 ! 195 209 IF(lwp) WRITE(numout,*) … … 399 413 ! ! ---------------------------------------------------- 400 414 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 401 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters415 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 402 416 DO jj = 2, jpjm1 403 417 DO ji = 2, jpim1 404 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji+1,jj) ) > &405 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &406 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) )&407 & 408 ll_tmp2 = ( ABS( sshn(ji+1,jj) -sshn(ji ,jj)) > 1.E-12 ).AND.( &409 & MAX( sshn(ji,jj) ,sshn(ji+1,jj) ) > &410 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )418 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 419 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 420 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 421 & > rn_wdmin1 + rn_wdmin2 422 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 423 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 424 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 411 425 412 426 IF(ll_tmp1) THEN … … 414 428 ELSE IF(ll_tmp2) THEN 415 429 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 417 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 430 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 431 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 432 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 433 418 434 ELSE 419 435 zcpx(ji,jj) = 0._wp 420 436 END IF 421 437 422 ll_tmp1 = MIN( sshn(ji,jj) ,sshn(ji,jj+1) ) > &423 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &424 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) )&425 & 426 ll_tmp2 = ( ABS( sshn(ji,jj) -sshn(ji,jj+1)) > 1.E-12 ).AND.( &427 & MAX( sshn(ji,jj) ,sshn(ji,jj+1) ) > &428 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )438 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 439 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 440 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 441 & > rn_wdmin1 + rn_wdmin2 442 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 443 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 444 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 429 445 430 446 IF(ll_tmp1) THEN … … 432 448 ELSE IF(ll_tmp2) THEN 433 449 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 434 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 435 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 450 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 451 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 452 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 453 436 454 ELSE 437 455 zcpy(ji,jj) = 0._wp … … 443 461 DO ji = 2, jpim1 444 462 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 445 & * r1_e1u(ji,jj) * zcpx(ji,jj) 463 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 446 464 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 447 & * r1_e2v(ji,jj) * zcpy(ji,jj) 465 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 466 448 467 END DO 449 468 END DO … … 468 487 END DO 469 488 ! 470 ! 489 ! ! Add bottom stress contribution from baroclinic velocities: 471 490 IF (ln_bt_fw) THEN 472 491 DO jj = 2, jpjm1 … … 490 509 ! 491 510 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 IF( ln_wd ) THEN493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 511 IF( ln_wd_il ) THEN 512 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) * wdrampu(ji,jj) 513 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) * wdrampv(ji,jj) 495 514 ELSE 496 515 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) … … 627 646 vn_adv(:,:) = 0._wp 628 647 ! ! ==================== ! 648 649 IF (ln_wd_dl) THEN 650 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 651 zvwdmask(:,:) = 0._wp ! 652 zuwdav2(:,:) = 0._wp 653 zvwdav2(:,:) = 0._wp 654 END IF 655 656 629 657 DO jn = 1, icycle ! sub-time-step loop ! 630 658 ! ! ==================== ! … … 654 682 ! Extrapolate Sea Level at step jit+0.5: 655 683 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 ! 684 685 ! set wetting & drying mask at tracer points for this barotropic sub-step 686 IF ( ln_wd_dl ) THEN 687 688 IF ( ln_wd_dl_rmp ) THEN 689 DO jj = 1, jpj 690 DO ji = 1, jpi ! vector opt. 691 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 692 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 693 ztwdmask(ji,jj) = 1._wp 694 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 695 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 696 ELSE 697 ztwdmask(ji,jj) = 0._wp 698 END IF 699 END DO 700 END DO 701 ELSE 702 DO jj = 1, jpj 703 DO ji = 1, jpi ! vector opt. 704 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 705 ztwdmask(ji,jj) = 1._wp 706 ELSE 707 ztwdmask(ji,jj) = 0._wp 708 END IF 709 END DO 710 END DO 711 END IF 712 713 END IF 714 715 657 716 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 658 717 DO ji = 2, fs_jpim1 ! Vector opt. … … 706 765 ENDIF 707 766 #endif 708 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 709 ! 767 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 768 769 IF ( ln_wd_dl ) THEN 770 771 772 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 773 774 DO jj = 1, jpjm1 775 DO ji = 1, jpim1 776 IF ( zwx(ji,jj) > 0.0 ) THEN 777 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 778 ELSE 779 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 780 END IF 781 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 782 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 783 784 IF ( zwy(ji,jj) > 0.0 ) THEN 785 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 786 ELSE 787 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 788 END IF 789 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 790 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 791 END DO 792 END DO 793 794 795 END IF 796 710 797 ! Sum over sub-time-steps to compute advective velocities 711 798 za2 = wgtbtp2(jn) 712 799 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 713 800 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 714 ! 801 802 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 803 IF ( ln_wd_dl_bc ) THEN 804 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 805 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 806 END IF 807 715 808 ! Set next sea level: 716 809 DO jj = 2, jpjm1 … … 766 859 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 767 860 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 768 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters861 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 769 862 DO jj = 2, jpjm1 770 863 DO ji = 2, jpim1 771 864 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 772 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &773 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &865 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 866 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 774 867 & > rn_wdmin1 + rn_wdmin2 775 868 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 776 869 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 777 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )870 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 778 871 779 872 IF(ll_tmp1) THEN … … 781 874 ELSE IF(ll_tmp2) THEN 782 875 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 783 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &876 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 784 877 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 785 878 ELSE … … 788 881 789 882 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 790 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &791 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &883 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 884 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 792 885 & > rn_wdmin1 + rn_wdmin2 793 886 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 794 887 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 795 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )888 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 796 889 797 890 IF(ll_tmp1) THEN … … 799 892 ELSE IF(ll_tmp2) THEN 800 893 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 801 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &894 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 802 895 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 803 896 ELSE … … 809 902 ! 810 903 ! Compute associated depths at U and V points: 811 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN 904 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 812 905 ! 813 906 DO jj = 2, jpjm1 … … 886 979 ! 887 980 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 981 !jth do implicitly instead 982 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 983 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 984 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 985 ENDIF 890 986 ! 891 987 ! Add top stresses: … … 895 991 ! Surface pressure trend: 896 992 897 IF( ln_wd ) THEN993 IF( ln_wd_il ) THEN 898 994 DO jj = 2, jpjm1 899 995 DO ji = 2, jpim1 … … 919 1015 ! 920 1016 ! Set next velocities: 921 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1017 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 922 1018 DO jj = 2, jpjm1 923 1019 DO ji = fs_2, fs_jpim1 ! vector opt. … … 933 1029 & + zv_frc(ji,jj) ) & 934 1030 & ) * ssvmask(ji,jj) 1031 1032 !jth implicit bottom friction: 1033 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1034 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 1035 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 1036 ENDIF 1037 935 1038 END DO 936 1039 END DO 937 1040 ! 938 ELSE 1041 ELSE !* Flux form 939 1042 DO jj = 2, jpjm1 940 1043 DO ji = fs_2, fs_jpim1 ! vector opt. 941 1044 942 IF( ln_wd ) THEN 943 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 944 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 945 ELSE 946 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 947 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 948 END IF 1045 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1046 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1047 949 1048 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 950 1049 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 964 1063 END DO 965 1064 ENDIF 966 ! 1065 1066 967 1067 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 968 IF( ln_wd ) THEN 969 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 970 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 971 ELSE 972 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 973 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 974 END IF 1068 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1069 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 975 1070 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 976 1071 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 1005 1100 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 1006 1101 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 1007 ELSE ! Sum transports 1008 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1009 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1010 ENDIF 1011 ! ! Sum sea level 1102 ELSE ! Sum transports 1103 IF ( .NOT.ln_wd_dl ) THEN 1104 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1105 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1106 ELSE 1107 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 1108 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 1109 END IF 1110 ENDIF 1111 ! ! Sum sea level 1012 1112 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 1113 1013 1114 ! ! ==================== ! 1014 1115 END DO ! end loop ! … … 1033 1134 vb2_b(:,:) = zwy(:,:) 1034 1135 ENDIF 1136 1137 1035 1138 ! 1036 1139 ! Update barotropic trend: … … 1062 1165 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1063 1166 ENDIF 1064 ! 1167 1168 1169 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1065 1170 DO jk = 1, jpkm1 1066 ! Correct velocities:1067 1171 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1068 1172 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1069 !1070 1173 END DO 1071 ! 1174 1175 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 1176 DO jk = 1, jpkm1 1177 un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk) 1178 vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk) 1179 END DO 1180 END IF 1181 1182 1072 1183 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1073 1184 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current … … 1097 1208 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1098 1209 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1210 IF( ln_wd_il ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1211 IF( ln_wd_dl ) CALL wrk_dealloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 1100 1212 ! 1101 1213 IF ( ln_diatmb ) THEN … … 1271 1383 rdtbt = rdt / REAL( nn_baro , wp ) 1272 1384 zcmax = zcmax * rdtbt 1273 1385 ! Print results 1274 1386 IF(lwp) WRITE(numout,*) 1275 1387 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
Note: See TracChangeset
for help on using the changeset viewer.
