Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 11 edited
-
dynadv.F90 (modified) (2 diffs)
-
dynhpg.F90 (modified) (18 diffs)
-
dynkeg.F90 (modified) (4 diffs)
-
dynldf.F90 (modified) (2 diffs)
-
dynnxt.F90 (modified) (3 diffs)
-
dynspg.F90 (modified) (3 diffs)
-
dynspg_ts.F90 (modified) (24 diffs)
-
dynvor.F90 (modified) (39 diffs)
-
dynzdf.F90 (modified) (1 diff)
-
sshwzv.F90 (modified) (5 diffs)
-
wet_dry.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r6140 r7646 106 106 WRITE(numout,*) 107 107 WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme' 108 WRITE(numout,*) '~~~~~~~~~~~ '109 WRITE(numout,*) ' 110 WRITE(numout,*) ' 111 WRITE(numout,*) ' 112 WRITE(numout,*) ' 113 WRITE(numout,*) ' 114 WRITE(numout,*) ' 108 WRITE(numout,*) '~~~~~~~~~~~~' 109 WRITE(numout,*) ' Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 110 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 111 WRITE(numout,*) ' = 0 standard scheme ; =1 Hollingsworth scheme nn_dynkeg = ', nn_dynkeg 112 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 113 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 114 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 115 115 ENDIF 116 116 … … 134 134 IF(lwp) THEN ! Print the choice 135 135 WRITE(numout,*) 136 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used'137 IF( nadv == 1 ) WRITE(numout,*) ' vector form : keg + zad_zts + vor is used'136 IF( nadv == 0 ) WRITE(numout,*) ' ===>> vector form : keg + zad + vor is used' 137 IF( nadv == 1 ) WRITE(numout,*) ' ===>> vector form : keg + zad_zts + vor is used' 138 138 IF( nadv == 0 .OR. nadv == 1 ) THEN 139 IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) ' with Centered standard keg scheme'140 IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) ' with Hollingsworth keg scheme'139 IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) ' with Centered standard keg scheme' 140 IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) ' with Hollingsworth keg scheme' 141 141 ENDIF 142 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used'143 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used'142 IF( nadv == 2 ) WRITE(numout,*) ' ===>> flux form : 2nd order scheme is used' 143 IF( nadv == 3 ) WRITE(numout,*) ' ===>> flux form : UBS scheme is used' 144 144 ENDIF 145 145 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6152 r7646 432 432 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 433 433 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 LOGICAL :: ll_tmp1, ll_tmp2 , ll_tmp3! local logical variables434 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 435 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter … … 438 438 ! 439 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )440 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 441 441 ! 442 442 IF( kt == nit000 ) THEN … … 451 451 ENDIF 452 452 ! 453 IF( ln_wd) THEN453 IF( ln_wd ) THEN 454 454 DO jj = 2, jpjm1 455 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 459 & rn_wdmin1 + rn_wdmin2 460 461 IF(ll_tmp1.AND.ll_tmp2) THEN 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 461 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 462 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 464 IF(ll_tmp1) THEN 462 465 zcpx(ji,jj) = 1.0_wp 463 wduflt(ji,jj) = 1.0_wp 464 ELSE IF(ll_tmp3) THEN 465 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 466 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 wduflt(ji,jj) = 1.0_wp 466 ELSE IF(ll_tmp2) THEN 467 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 469 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 469 470 ELSE 470 471 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp472 472 END IF 473 473 474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 477 & rn_wdmin1 + rn_wdmin2 478 479 IF(ll_tmp1.AND.ll_tmp2) THEN 474 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 476 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 477 & > rn_wdmin1 + rn_wdmin2 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 479 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 480 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 482 IF(ll_tmp1) THEN 480 483 zcpy(ji,jj) = 1.0_wp 481 wdvflt(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp3) THEN 483 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 wdvflt(ji,jj) = 1.0_wp 484 ELSE IF(ll_tmp2) THEN 485 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 487 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 487 488 ELSE 488 489 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp490 490 END IF 491 491 END DO 492 492 END DO 493 493 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 494 ENDIF 495 494 END IF 496 495 497 496 ! Surface value … … 510 509 511 510 512 IF( ln_wd) THEN511 IF( ln_wd ) THEN 513 512 514 513 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 541 540 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 542 541 543 IF( ln_wd) THEN542 IF( ln_wd ) THEN 544 543 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 545 544 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 556 555 ! 557 556 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 558 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )557 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 559 558 ! 560 559 END SUBROUTINE hpg_sco … … 701 700 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 702 701 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 703 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )704 ! 705 ! 706 IF( ln_wd) THEN702 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 703 ! 704 ! 705 IF( ln_wd ) THEN 707 706 DO jj = 2, jpjm1 708 707 DO ji = 2, jpim1 709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 713 & rn_wdmin1 + rn_wdmin2 708 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 709 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 710 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 713 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 714 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 714 715 715 716 IF(ll_tmp1) THEN 716 717 zcpx(ji,jj) = 1.0_wp 717 718 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here719 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&720 & (sshn(ji+1,jj) - sshn(ji,jj)))719 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 721 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 721 722 ELSE 722 723 zcpx(ji,jj) = 0._wp 723 724 END IF 724 725 725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 729 & rn_wdmin1 + rn_wdmin2 726 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 727 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 728 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 729 & > rn_wdmin1 + rn_wdmin2 730 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 731 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 732 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 730 733 731 734 IF(ll_tmp1) THEN 732 735 zcpy(ji,jj) = 1.0_wp 733 736 ELSE IF(ll_tmp2) THEN 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here735 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&736 & (sshn(ji,jj+1) - sshn(ji,jj)))737 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 739 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 737 740 ELSE 738 741 zcpy(ji,jj) = 0._wp … … 741 744 END DO 742 745 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 743 ENDIF 744 746 END IF 745 747 746 748 IF( kt == nit000 ) THEN … … 913 915 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 914 916 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 915 IF( ln_wd) THEN917 IF( ln_wd ) THEN 916 918 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 917 919 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) … … 936 938 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 937 939 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 938 IF( ln_wd) THEN940 IF( ln_wd ) THEN 939 941 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 940 942 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 950 952 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 951 953 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 952 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )954 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 953 955 ! 954 956 END SUBROUTINE hpg_djc … … 987 989 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 988 990 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 989 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )991 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 990 992 ! 991 993 IF( kt == nit000 ) THEN … … 1000 1002 IF( ln_linssh ) znad = 0._wp 1001 1003 1002 IF( ln_wd) THEN1004 IF( ln_wd ) THEN 1003 1005 DO jj = 2, jpjm1 1004 1006 DO ji = 2, jpim1 1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1009 & rn_wdmin1 + rn_wdmin2 1007 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 1009 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 1010 & > rn_wdmin1 + rn_wdmin2 1011 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 1012 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1013 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1010 1014 1011 1015 IF(ll_tmp1) THEN 1012 1016 zcpx(ji,jj) = 1.0_wp 1013 1017 ELSE IF(ll_tmp2) THEN 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here1015 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&1016 & (sshn(ji+1,jj) - sshn(ji,jj)))1018 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 1020 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1017 1021 ELSE 1018 1022 zcpx(ji,jj) = 0._wp 1019 1023 END IF 1020 1024 1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1025 & rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1.OR.ll_tmp2) THEN 1025 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1026 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 1027 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 1028 & > rn_wdmin1 + rn_wdmin2 1029 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 1030 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1031 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 1033 IF(ll_tmp1) THEN 1028 1034 zcpy(ji,jj) = 1.0_wp 1029 1035 ELSE IF(ll_tmp2) THEN 1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here1031 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&1032 & (sshn(ji,jj+1) - sshn(ji,jj)))1036 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1037 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 1038 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1033 1039 ELSE 1034 1040 zcpy(ji,jj) = 0._wp … … 1037 1043 END DO 1038 1044 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 END IF1045 END IF 1040 1046 1041 1047 ! Clean 3-D work arrays … … 1046 1052 DO jj = 1, jpj 1047 1053 DO ji = 1, jpi 1048 jk = mb athy(ji,jj)1054 jk = mbkt(ji,jj)+1 1049 1055 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 1056 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) … … 1221 1227 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1222 1228 ENDIF 1223 IF( ln_wd) THEN1229 IF( ln_wd ) THEN 1224 1230 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 1231 zdpdx2 = zdpdx2 * zcpx(ji,jj) … … 1280 1286 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1281 1287 ENDIF 1282 IF( ln_wd) THEN1288 IF( ln_wd ) THEN 1283 1289 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 1290 zdpdy2 = zdpdy2 * zcpy(ji,jj) … … 1295 1301 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1296 1302 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )1303 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1298 1304 ! 1299 1305 END SUBROUTINE hpg_prj -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5328 r7646 24 24 USE wrk_nemo ! Memory Allocation 25 25 USE timing ! Timing 26 USE bdy_oce ! ocean open boundary conditions 26 27 27 28 IMPLICIT NONE … … 78 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 79 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 81 INTEGER :: jb ! dummy loop indices 82 INTEGER :: ii, ij, igrd, ib_bdy ! local integers 83 INTEGER :: fu, fv 80 84 !!---------------------------------------------------------------------- 81 85 ! … … 98 102 zhke(:,:,jpk) = 0._wp 99 103 104 IF (ln_bdy) THEN 105 ! Maria Luneva & Fred Wobus: July-2016 106 ! compensate for lack of turbulent kinetic energy on liquid bdy points 107 DO ib_bdy = 1, nb_bdy 108 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 109 igrd = 2 ! Copying normal velocity into points outside bdy 110 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 111 DO jk = 1, jpkm1 112 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 113 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 114 fu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 115 un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 116 END DO 117 END DO 118 ! 119 igrd = 3 ! Copying normal velocity into points outside bdy 120 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 121 DO jk = 1, jpkm1 122 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 123 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 124 fv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 125 vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 126 END DO 127 END DO 128 ENDIF 129 ENDDO 130 ENDIF 131 100 132 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 133 ! … … 133 165 ! 134 166 END SELECT 167 168 IF (ln_bdy) THEN 169 ! restore velocity masks at points outside boundary 170 un(:,:,:) = un(:,:,:) * umask(:,:,:) 171 vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 172 ENDIF 173 174 135 175 ! 136 176 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r6140 r7646 110 110 WRITE(numout,*) 111 111 WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics' 112 WRITE(numout,*) '~~~~~~~~~~~ '113 WRITE(numout,*) ' 114 WRITE(numout,*) ' 115 WRITE(numout,*) ' 116 WRITE(numout,*) ' 117 WRITE(numout,*) ' 118 WRITE(numout,*) ' 112 WRITE(numout,*) '~~~~~~~~~~~~' 113 WRITE(numout,*) ' Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 114 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 115 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 116 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 117 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 118 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 119 119 ENDIF 120 120 ! ! use of lateral operator or not … … 180 180 IF(lwp) THEN 181 181 WRITE(numout,*) 182 IF( nldf == np_no_ldf ) WRITE(numout,*) ' 183 IF( nldf == np_lap ) WRITE(numout,*) ' 184 IF( nldf == np_lap_i ) WRITE(numout,*) ' 185 IF( nldf == np_blp ) WRITE(numout,*) ' 182 IF( nldf == np_no_ldf ) WRITE(numout,*) ' ===>> NO lateral viscosity' 183 IF( nldf == np_lap ) WRITE(numout,*) ' ===>> iso-level laplacian operator' 184 IF( nldf == np_lap_i ) WRITE(numout,*) ' ===>> rotated laplacian operator with iso-level background' 185 IF( nldf == np_blp ) WRITE(numout,*) ' ===>> iso-level bi-laplacian operator' 186 186 ENDIF 187 187 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r6140 r7646 32 32 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 33 33 USE domvvl ! variable volume 34 USE bdy_oce ! ocean open boundary conditions34 USE bdy_oce , ONLY: ln_bdy 35 35 USE bdydta ! ocean open boundary conditions 36 36 USE bdydyn ! ocean open boundary conditions … … 77 77 !! * Apply lateral boundary conditions on after velocity 78 78 !! at the local domain boundaries through lbc_lnk call, 79 !! at the one-way open boundaries (l k_bdy=T),79 !! at the one-way open boundaries (ln_bdy=T), 80 80 !! at the AGRIF zoom boundaries (lk_agrif=T) 81 81 !! … … 147 147 CALL lbc_lnk( va, 'V', -1. ) 148 148 ! 149 # if defined key_bdy150 149 ! !* BDY open boundaries 151 IF( l k_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )152 IF( l k_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )150 IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 151 IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 153 152 154 153 !!$ Do we need a call to bdy_vol here?? 155 !156 # endif157 154 ! 158 155 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r6981 r7646 88 88 ! 89 89 IF( ln_apr_dyn & ! atmos. pressure 90 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. l k_tide) ) & ! tide potential (no time slitting)90 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) ) & ! tide potential (no time slitting) 91 91 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 92 92 ! … … 111 111 ! 112 112 ! !== tide potential forcing term ==! 113 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. l k_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case113 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 114 114 ! 115 115 CALL upd_tide( kt ) ! update tide potential … … 216 216 IF(lwp) THEN 217 217 WRITE(numout,*) 218 IF( nspg == np_EXP ) WRITE(numout,*) ' explicit free surface'219 IF( nspg == np_TS ) WRITE(numout,*) ' free surface with time splitting scheme'220 IF( nspg == np_NO ) WRITE(numout,*) ' No surface surface pressure gradient trend in momentum Eqs.'218 IF( nspg == np_EXP ) WRITE(numout,*) ' ===>> explicit free surface' 219 IF( nspg == np_TS ) WRITE(numout,*) ' ===>> free surface with time splitting scheme' 220 IF( nspg == np_NO ) WRITE(numout,*) ' ===>> No surface surface pressure gradient trend in momentum Eqs.' 221 221 ENDIF 222 222 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r6152 r7646 15 15 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 16 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 17 !! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence 17 18 !!--------------------------------------------------------------------- 18 19 … … 33 34 USE dynvor ! vorticity term 34 35 USE wet_dry ! wetting/drying flux limter 35 USE bdy_ par ! for lk_bdy36 USE bdy_oce ! open boundary 36 37 USE bdytides ! open boundary condition data 37 38 USE bdydyn2d ! open boundary conditions on barotropic variables 38 39 USE sbctide ! tides 39 40 USE updtide ! tide potential 41 USE sbcwave ! surface wave 40 42 ! 41 43 USE in_out_manager ! I/O manager … … 69 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 70 72 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff /h at F points73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff_f/h at F points 72 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 73 75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) … … 156 158 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 159 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef.159 160 !!---------------------------------------------------------------------- 160 161 ! … … 168 169 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 169 170 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 171 172 ! 172 173 zmdi=1.e+20 ! missing data indicator for masking … … 220 221 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 221 222 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 222 SELECT CASE( nn_een_e3f ) !* ff /e3 at F-point223 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 223 224 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 224 225 DO jj = 1, jpjm1 … … 226 227 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 227 228 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 228 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff (ji,jj) / zwz(ji,jj)229 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 229 230 END DO 230 231 END DO … … 236 237 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 237 238 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 238 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff (ji,jj) / zwz(ji,jj)239 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 239 240 END DO 240 241 END DO … … 255 256 zwz(:,:) = 0._wp 256 257 zhf(:,:) = 0._wp 257 IF ( .not. ln_sco ) THEN 258 259 !!gm agree the JC comment : this should be done in a much clear way 260 261 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 262 ! Set it to zero for the time being 263 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 264 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 265 ! ENDIF 266 ! zhf(:,:) = gdepw_0(:,:,jk+1) 267 ELSE 268 zhf(:,:) = hbatf(:,:) 269 END IF 270 271 DO jj = 1, jpjm1 272 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 273 END DO 258 259 !!gm assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed 260 !!gm A priori a better value should be something like : 261 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 262 !!gm divided by the sum of the corresponding mask 263 !!gm 264 !! 265 IF ( .not. ln_sco ) THEN 266 267 !!gm agree the JC comment : this should be done in a much clear way 268 269 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 270 ! Set it to zero for the time being 271 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 272 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 273 ! ENDIF 274 ! zhf(:,:) = gdepw_0(:,:,jk+1) 275 ELSE 276 !zhf(:,:) = hbatf(:,:) 277 DO jj = 1, jpjm1 278 DO ji = 1, jpim1 279 zhf(ji,jj) = MAX( 0._wp, & 280 & ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + & 281 & ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + & 282 & ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + & 283 & ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / & 284 & ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +& 285 & tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +& 286 & rsmall ) ) 287 END DO 288 END DO 289 END IF 290 291 DO jj = 1, jpjm1 292 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 293 END DO 294 !!gm end 274 295 275 296 DO jk = 1, jpkm1 … … 285 306 END DO 286 307 END DO 287 zwz(:,:) = ff (:,:) * zwz(:,:)308 zwz(:,:) = ff_f(:,:) * zwz(:,:) 288 309 ENDIF 289 310 ENDIF … … 324 345 END DO 325 346 END DO 347 348 !!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 349 !!gm Is it correct to do so ? I think so... 350 351 326 352 ! !* barotropic Coriolis trends (vorticity scheme dependent) 327 353 ! ! -------------------------------------------------------- … … 374 400 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 375 401 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp 377 wdvflt1(:,:) = 1.0_wp 378 DO jj = 2, jpjm1 379 DO ji = 2, jpim1 380 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 381 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 382 & > rn_wdmin1 + rn_wdmin2 383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 384 & + rn_wdmin1 + rn_wdmin2 402 DO jj = 2, jpjm1 403 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 & > rn_wdmin1 + rn_wdmin2 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 ) 411 385 412 IF(ll_tmp1) THEN 386 zcpx(ji,jj) 387 ELSE IF(ll_tmp2) THEN388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happenhere389 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) &390 & /(sshn(ji+1,jj) - sshn(ji,jj)))413 zcpx(ji,jj) = 1.0_wp 414 ELSE IF(ll_tmp2) THEN 415 ! 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)) ) 391 418 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 419 zcpx(ji,jj) = 0._wp 394 420 END IF 395 396 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 397 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 398 & > rn_wdmin1 + rn_wdmin2 399 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 400 & + rn_wdmin1 + rn_wdmin2 421 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 & > rn_wdmin1 + rn_wdmin2 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 ) 429 401 430 IF(ll_tmp1) THEN 402 zcpy(ji,jj)= 1.0_wp403 ELSE IF(ll_tmp2) THEN404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happenhere405 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) &406 & /(sshn(ji,jj+1) - sshn(ji,jj)))431 zcpy(ji,jj) = 1.0_wp 432 ELSE IF(ll_tmp2) THEN 433 ! 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 )) ) 407 436 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 437 zcpy(ji,jj) = 0._wp 438 END IF 439 END DO 413 440 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 441 417 442 DO jj = 2, jpjm1 418 443 DO ji = 2, jpim1 419 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &420 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj)421 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &422 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj)444 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 445 & * r1_e1u(ji,jj) * zcpx(ji,jj) 446 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 447 & * r1_e2v(ji,jj) * zcpy(ji,jj) 423 448 END DO 424 449 END DO … … 474 499 ! 475 500 ! ! Add top stress contribution from baroclinic velocities: 476 IF (ln_bt_fw) THEN501 IF( ln_bt_fw ) THEN 477 502 DO jj = 2, jpjm1 478 503 DO ji = fs_2, fs_jpim1 ! vector opt. … … 538 563 & + fwfisf(:,:) + fwfisf_b(:,:) ) 539 564 ENDIF 565 ! 566 IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary 567 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 568 ENDIF 569 ! 540 570 #if defined key_asminc 541 571 ! ! Include the IAU weighted SSH increment … … 567 597 ENDIF 568 598 569 IF( ln_wd ) THEN !preserve the positivity of water depth570 !ssh[b,n,a] should have already been processed for this571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))573 ENDIF574 599 ! 575 600 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 608 633 ! Update only tidal forcing at open boundaries 609 634 #if defined key_tide 610 IF( l k_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )611 IF( ln_tide_pot .AND. l k_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )635 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 636 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 612 637 #endif 613 638 ! … … 646 671 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 647 672 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 648 IF( ln_wd ) THEN649 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)650 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)651 END IF652 673 ELSE 653 674 zhup2_e (:,:) = hu_n(:,:) … … 702 723 END DO 703 724 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 704 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))725 705 726 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 706 727 707 #if defined key_bdy708 728 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 709 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 710 #endif 729 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) 711 730 #if defined key_agrif 712 731 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) … … 750 769 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 751 770 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 752 wduflt1(:,:) = 1._wp753 wdvflt1(:,:) = 1._wp754 771 DO jj = 2, jpjm1 755 DO ji = 2, jpim1 756 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 757 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 758 & > rn_wdmin1 + rn_wdmin2 759 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 760 & + rn_wdmin1 + rn_wdmin2 761 IF(ll_tmp1) THEN 762 zcpx(ji,jj) = 1._wp 763 ELSE IF(ll_tmp2) THEN 764 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 765 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 766 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 767 ELSE 768 zcpx(ji,jj) = 0._wp 769 wduflt1(ji,jj) = 0._wp 770 END IF 771 772 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 773 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 776 & + rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpy(ji,jj) = 1._wp 779 ELSE IF(ll_tmp2) THEN 780 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 781 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 782 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 783 ELSE 784 zcpy(ji,jj) = 0._wp 785 wdvflt1(ji,jj) = 0._wp 786 END IF 772 DO ji = 2, jpim1 773 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 774 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 775 & MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) & 776 & > rn_wdmin1 + rn_wdmin2 777 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 778 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 779 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 780 781 IF(ll_tmp1) THEN 782 zcpx(ji,jj) = 1.0_wp 783 ELSE IF(ll_tmp2) THEN 784 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 785 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 786 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 787 ELSE 788 zcpx(ji,jj) = 0._wp 789 END IF 790 791 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 792 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 793 & MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) & 794 & > rn_wdmin1 + rn_wdmin2 795 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 796 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 797 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 798 799 IF(ll_tmp1) THEN 800 zcpy(ji,jj) = 1.0_wp 801 ELSE IF(ll_tmp2) THEN 802 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 803 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 804 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 805 ELSE 806 zcpy(ji,jj) = 0._wp 807 END IF 787 808 END DO 788 END DO 789 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 790 ENDIF 809 END DO 810 END IF 791 811 ! 792 812 ! Compute associated depths at U and V points: … … 806 826 END DO 807 827 808 IF( ln_wd ) THEN809 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )810 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )811 END IF812 813 828 ENDIF 814 829 ! … … 861 876 ! 862 877 ! Add tidal astronomical forcing if defined 863 IF ( l k_tide.AND.ln_tide_pot ) THEN878 IF ( ln_tide .AND. ln_tide_pot ) THEN 864 879 DO jj = 2, jpjm1 865 880 DO ji = fs_2, fs_jpim1 ! vector opt. … … 967 982 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 968 983 ! 969 #if defined key_bdy970 984 ! ! open boundaries 971 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 972 #endif 985 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 973 986 #if defined key_agrif 974 987 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif … … 1086 1099 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1087 1100 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)1101 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1089 1102 ! 1090 1103 IF ( ln_diatmb ) THEN … … 1248 1261 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1249 1262 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1250 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) )1263 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 1251 1264 END DO 1252 1265 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r6140 r7646 17 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 18 18 !! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory 19 !! - ! 2016-12 (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 19 20 !!---------------------------------------------------------------------- 20 21 … … 32 33 USE trd_oce ! trends: ocean variables 33 34 USE trddyn ! trend manager: dynamics 35 USE sbcwave ! Surface Waves (add Stokes-Coriolis force) 36 USE sbc_oce , ONLY : ln_stcor ! use Stoke-Coriolis force 34 37 ! 35 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 77 80 # include "vectopt_loop_substitute.h90" 78 81 !!---------------------------------------------------------------------- 79 !! NEMO/OPA 3.7 , NEMO Consortium (201 4)82 !! NEMO/OPA 3.7 , NEMO Consortium (2016) 80 83 !! $Id$ 81 84 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 108 111 ztrdu(:,:,:) = ua(:,:,:) 109 112 ztrdv(:,:,:) = va(:,:,:) 110 CALL vor_ene( kt, nrvm, u a, va )! relative vorticity or metric trend113 CALL vor_ene( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend 111 114 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 112 115 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 114 117 ztrdu(:,:,:) = ua(:,:,:) 115 118 ztrdv(:,:,:) = va(:,:,:) 116 CALL vor_ene( kt, ncor, u a, va )! planetary vorticity trend119 CALL vor_ene( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend 117 120 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 118 121 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 119 122 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 120 ELSE 121 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend 123 ELSE ! total vorticity trend 124 CALL vor_ene( kt, ntot, un , vn , ua, va ) ! total vorticity trend 125 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 122 126 ENDIF 123 127 ! … … 126 130 ztrdu(:,:,:) = ua(:,:,:) 127 131 ztrdv(:,:,:) = va(:,:,:) 128 CALL vor_ens( kt, nrvm, u a, va )! relative vorticity or metric trend132 CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend 129 133 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 130 134 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 132 136 ztrdu(:,:,:) = ua(:,:,:) 133 137 ztrdv(:,:,:) = va(:,:,:) 134 CALL vor_ens( kt, ncor, u a, va )! planetary vorticity trend138 CALL vor_ens( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend 135 139 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 136 140 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 141 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 ELSE 139 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend 142 ELSE ! total vorticity trend 143 CALL vor_ens( kt, ntot, un , vn , ua, va ) ! total vorticity trend 144 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 140 145 ENDIF 141 146 ! … … 144 149 ztrdu(:,:,:) = ua(:,:,:) 145 150 ztrdv(:,:,:) = va(:,:,:) 146 CALL vor_ens( kt, nrvm, u a, va )! relative vorticity or metric trend (ens)151 CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend (ens) 147 152 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 153 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 150 155 ztrdu(:,:,:) = ua(:,:,:) 151 156 ztrdv(:,:,:) = va(:,:,:) 152 CALL vor_ene( kt, ncor, u a, va )! planetary vorticity trend (ene)157 CALL vor_ene( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend (ene) 153 158 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 154 159 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 155 160 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 ELSE 157 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 158 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 161 ELSE ! total vorticity trend 162 CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend (ens) 163 CALL vor_ene( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend (ene) 164 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 159 165 ENDIF 160 166 ! … … 163 169 ztrdu(:,:,:) = ua(:,:,:) 164 170 ztrdv(:,:,:) = va(:,:,:) 165 CALL vor_een( kt, nrvm, u a, va )! relative vorticity or metric trend171 CALL vor_een( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend 166 172 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 167 173 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 169 175 ztrdu(:,:,:) = ua(:,:,:) 170 176 ztrdv(:,:,:) = va(:,:,:) 171 CALL vor_een( kt, ncor, u a, va )! planetary vorticity trend177 CALL vor_een( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend 172 178 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 173 179 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 174 180 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 175 ELSE 176 CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend 181 ELSE ! total vorticity trend 182 CALL vor_een( kt, ntot, un , vn , ua, va ) ! total vorticity trend 183 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 177 184 ENDIF 178 185 ! … … 190 197 191 198 192 SUBROUTINE vor_ene( kt, kvor, pu a, pva )199 SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 193 200 !!---------------------------------------------------------------------- 194 201 !! *** ROUTINE vor_ene *** … … 210 217 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 211 218 !!---------------------------------------------------------------------- 212 INTEGER , INTENT(in ) :: kt ! ocean time-step index213 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;214 ! ! =nrvm (relative vorticity or metric)215 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu a ! total u-trend216 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: p va ! total v-trend219 INTEGER , INTENT(in ) :: kt ! ocean time-step index 220 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 221 ! ! =nrvm (relative vorticity or metric) 222 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! now velocities 223 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua, pva ! total v-trend 217 224 ! 218 225 INTEGER :: ji, jj, jk ! dummy loop indices … … 223 230 IF( nn_timing == 1 ) CALL timing_start('vor_ene') 224 231 ! 225 CALL wrk_alloc( jpi, jpj,zwx, zwy, zwz )232 CALL wrk_alloc( jpi,jpj, zwx, zwy, zwz ) 226 233 ! 227 234 IF( kt == nit000 ) THEN … … 237 244 SELECT CASE( kvor ) !== vorticity considered ==! 238 245 CASE ( np_COR ) !* Coriolis (planetary vorticity) 239 zwz(:,:) = ff (:,:)246 zwz(:,:) = ff_f(:,:) 240 247 CASE ( np_RVO ) !* relative vorticity 241 248 DO jj = 1, jpjm1 242 249 DO ji = 1, fs_jpim1 ! vector opt. 243 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &244 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) * r1_e1e2f(ji,jj)250 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 251 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 245 252 END DO 246 253 END DO … … 248 255 DO jj = 1, jpjm1 249 256 DO ji = 1, fs_jpim1 ! vector opt. 250 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &251 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &257 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 258 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 252 259 & * 0.5 * r1_e1e2f(ji,jj) 253 260 END DO … … 256 263 DO jj = 1, jpjm1 257 264 DO ji = 1, fs_jpim1 ! vector opt. 258 zwz(ji,jj) = ff (ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &259 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) &265 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 266 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 260 267 & * r1_e1e2f(ji,jj) 261 268 END DO … … 264 271 DO jj = 1, jpjm1 265 272 DO ji = 1, fs_jpim1 ! vector opt. 266 zwz(ji,jj) = ff (ji,jj) &267 & + ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &268 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &273 zwz(ji,jj) = ff_f(ji,jj) & 274 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 275 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 269 276 & * 0.5 * r1_e1e2f(ji,jj) 270 277 END DO … … 284 291 IF( ln_sco ) THEN 285 292 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 286 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)287 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)293 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 294 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 288 295 ELSE 289 zwx(:,:) = e2u(:,:) * un(:,:,jk)290 zwy(:,:) = e1v(:,:) * vn(:,:,jk)296 zwx(:,:) = e2u(:,:) * pun(:,:,jk) 297 zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 291 298 ENDIF 292 299 ! !== compute and add the vorticity term trend =! … … 311 318 312 319 313 SUBROUTINE vor_ens( kt, kvor, pu a, pva )320 SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 314 321 !!---------------------------------------------------------------------- 315 322 !! *** ROUTINE vor_ens *** … … 331 338 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 332 339 !!---------------------------------------------------------------------- 333 INTEGER , INTENT(in ) :: kt ! ocean time-step index334 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;335 ! ! =nrvm (relative vorticity or metric)336 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu a ! total u-trend337 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: p va ! total v-trend340 INTEGER , INTENT(in ) :: kt ! ocean time-step index 341 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 342 ! ! =nrvm (relative vorticity or metric) 343 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! now velocities 344 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua, pva ! total v-trend 338 345 ! 339 346 INTEGER :: ji, jj, jk ! dummy loop indices … … 344 351 IF( nn_timing == 1 ) CALL timing_start('vor_ens') 345 352 ! 346 CALL wrk_alloc( jpi, jpj,zwx, zwy, zwz )353 CALL wrk_alloc( jpi,jpj, zwx, zwy, zwz ) 347 354 ! 348 355 IF( kt == nit000 ) THEN … … 357 364 SELECT CASE( kvor ) !== vorticity considered ==! 358 365 CASE ( np_COR ) !* Coriolis (planetary vorticity) 359 zwz(:,:) = ff (:,:)366 zwz(:,:) = ff_f(:,:) 360 367 CASE ( np_RVO ) !* relative vorticity 361 368 DO jj = 1, jpjm1 362 369 DO ji = 1, fs_jpim1 ! vector opt. 363 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &364 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) * r1_e1e2f(ji,jj)370 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 371 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 365 372 END DO 366 373 END DO … … 368 375 DO jj = 1, jpjm1 369 376 DO ji = 1, fs_jpim1 ! vector opt. 370 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &371 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &377 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 378 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 372 379 & * 0.5 * r1_e1e2f(ji,jj) 373 380 END DO … … 376 383 DO jj = 1, jpjm1 377 384 DO ji = 1, fs_jpim1 ! vector opt. 378 zwz(ji,jj) = ff (ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &379 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) &385 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 386 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 380 387 & * r1_e1e2f(ji,jj) 381 388 END DO … … 384 391 DO jj = 1, jpjm1 385 392 DO ji = 1, fs_jpim1 ! vector opt. 386 zwz(ji,jj) = ff (ji,jj) &387 & + ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &388 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &393 zwz(ji,jj) = ff_f(ji,jj) & 394 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 395 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 389 396 & * 0.5 * r1_e1e2f(ji,jj) 390 397 END DO … … 404 411 IF( ln_sco ) THEN !== horizontal fluxes ==! 405 412 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 406 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)407 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)413 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 414 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 408 415 ELSE 409 zwx(:,:) = e2u(:,:) * un(:,:,jk)410 zwy(:,:) = e1v(:,:) * vn(:,:,jk)416 zwx(:,:) = e2u(:,:) * pun(:,:,jk) 417 zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 411 418 ENDIF 412 419 ! !== compute and add the vorticity term trend =! … … 431 438 432 439 433 SUBROUTINE vor_een( kt, kvor, pu a, pva )440 SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 434 441 !!---------------------------------------------------------------------- 435 442 !! *** ROUTINE vor_een *** … … 448 455 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 449 456 !!---------------------------------------------------------------------- 450 INTEGER , INTENT(in ) :: kt ! ocean time-step index 451 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 452 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 453 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 457 INTEGER , INTENT(in ) :: kt ! ocean time-step index 458 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 459 ! ! =nrvm (relative vorticity or metric) 460 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! now velocities 461 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua, pva ! total v-trend 454 462 ! 455 463 INTEGER :: ji, jj, jk ! dummy loop indices … … 506 514 DO jj = 1, jpjm1 507 515 DO ji = 1, fs_jpim1 ! vector opt. 508 zwz(ji,jj) = ff (ji,jj) * z1_e3f(ji,jj)516 zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj) 509 517 END DO 510 518 END DO … … 512 520 DO jj = 1, jpjm1 513 521 DO ji = 1, fs_jpim1 ! vector opt. 514 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &515 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) &522 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 523 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 516 524 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 517 525 END DO … … 520 528 DO jj = 1, jpjm1 521 529 DO ji = 1, fs_jpim1 ! vector opt. 522 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &523 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &530 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 531 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 524 532 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 525 533 END DO … … 528 536 DO jj = 1, jpjm1 529 537 DO ji = 1, fs_jpim1 ! vector opt. 530 zwz(ji,jj) = ( ff (ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) *vn(ji,jj,jk) &531 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) *un(ji,jj,jk) ) &538 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 539 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 532 540 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 533 541 END DO … … 536 544 DO jj = 1, jpjm1 537 545 DO ji = 1, fs_jpim1 ! vector opt. 538 zwz(ji,jj) = ( ff (ji,jj) &539 & + ( ( vn(ji+1,jj ,jk) +vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &540 & - ( un(ji ,jj+1,jk) +un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &546 zwz(ji,jj) = ( ff_f(ji,jj) & 547 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 548 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 541 549 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 542 550 END DO … … 557 565 ! 558 566 ! !== horizontal fluxes ==! 559 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)560 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)567 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 568 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 561 569 562 570 ! !== compute and add the vorticity term trend =! … … 626 634 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 627 635 WRITE(numout,*) '~~~~~~~~~~~~' 628 WRITE(numout,*) ' 629 WRITE(numout,*) ' 630 WRITE(numout,*) ' 631 WRITE(numout,*) ' 632 WRITE(numout,*) ' 633 WRITE(numout,*) ' 634 WRITE(numout,*) ' 636 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 637 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 638 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 639 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 640 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 641 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 642 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 635 643 ENDIF 636 644 … … 639 647 ! at angles with three ocean points and one land point 640 648 IF(lwp) WRITE(numout,*) 641 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T)ln_vorlat = ', ln_vorlat649 IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 642 650 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 643 651 DO jk = 1, jpk … … 666 674 ncor = np_COR 667 675 IF( ln_dynadv_vec ) THEN 668 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'676 IF(lwp) WRITE(numout,*) ' ===>> Vector form advection : vorticity = Coriolis + relative vorticity' 669 677 nrvm = np_RVO ! relative vorticity 670 678 ntot = np_CRV ! relative + planetary vorticity 671 679 ELSE 672 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'680 IF(lwp) WRITE(numout,*) ' ===>> Flux form advection : vorticity = Coriolis + metric term' 673 681 nrvm = np_MET ! metric term 674 682 ntot = np_CME ! Coriolis + metric term … … 677 685 IF(lwp) THEN ! Print the choice 678 686 WRITE(numout,*) 679 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>>energy conserving scheme'680 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>>enstrophy conserving scheme'681 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>>mixed enstrophy/energy conserving scheme'682 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>>energy and enstrophy conserving scheme'687 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' ===>> energy conserving scheme' 688 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' ===>> enstrophy conserving scheme' 689 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' ===>> mixed enstrophy/energy conserving scheme' 690 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' ===>> energy and enstrophy conserving scheme' 683 691 ENDIF 684 692 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r6140 r7646 119 119 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 120 120 WRITE(numout,*) '~~~~~~~~~~~' 121 IF( nzdf == 0 ) WRITE(numout,*) ' 122 IF( nzdf == 1 ) WRITE(numout,*) ' 121 IF( nzdf == 0 ) WRITE(numout,*) ' ===>> Explicit time-splitting scheme' 122 IF( nzdf == 1 ) WRITE(numout,*) ' ===>> Implicit (euler backward) scheme' 123 123 ENDIF 124 124 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r6152 r7646 22 22 USE divhor ! horizontal divergence 23 23 USE phycst ! physical constants 24 USE bdy_oce ! 25 USE bdy_par ! 24 USE bdy_oce , ONLY: ln_bdy, bdytmask 26 25 USE bdydyn2d ! bdy_ssh routine 27 26 #if defined key_agrif … … 88 87 ENDIF 89 88 ! 90 CALL div_hor( kt ) ! Horizontal divergence 91 ! 92 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 89 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 93 90 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 91 zcoef = 0.5_wp * r1_rau0 94 92 95 93 ! !------------------------------! 96 94 ! ! After Sea Surface Height ! 97 95 ! !------------------------------! 96 IF(ln_wd) THEN 97 CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 98 ENDIF 99 100 CALL div_hor( kt ) ! Horizontal divergence 101 ! 98 102 zhdiv(:,:) = 0._wp 99 103 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports … … 104 108 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 105 109 ! 106 zcoef = 0.5_wp * r1_rau0107 108 IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)109 110 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 111 … … 116 116 CALL agrif_ssh( kt ) 117 117 # endif 118 # if defined key_bdy 119 IF( lk_bdy ) THEN 118 IF( ln_bdy ) THEN 120 119 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 121 120 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 122 121 ENDIF 123 # endif124 122 ENDIF 125 123 … … 211 209 ENDIF 212 210 213 #if defined key_bdy 214 IF( lk_bdy ) THEN 211 IF( ln_bdy ) THEN 215 212 DO jk = 1, jpkm1 216 213 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 217 214 END DO 218 215 ENDIF 219 #endif220 216 ! 221 217 IF( nn_timing == 1 ) CALL timing_stop('wzv') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r6152 r7646 1 2 1 MODULE wet_dry 3 2 !!============================================================================== … … 7 6 !! only effects if wetting/drying is on (ln_wd == .true.) 8 7 !!============================================================================== 9 !! History : 10 !! NEMO 3.6 ! 2014-09 ((H.Liu) Original code 8 !! History : 3.6 ! 2014-09 ((H.Liu) Original code 11 9 !! ! will add the runoff and periodic BC case later 12 10 !!---------------------------------------------------------------------- … … 33 31 !! --------------------------------------------------------------------- 34 32 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter36 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths 35 ! (can include negative depths) 37 36 38 37 LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off) … … 77 76 WRITE(numout,*) 78 77 WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read' 79 WRITE(numout,*) '~~~~~~~ 78 WRITE(numout,*) '~~~~~~~~' 80 79 WRITE(numout,*) ' Namelist namwad' 81 80 WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd … … 84 83 WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld 85 84 WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit 86 87 85 ENDIF 86 ! 88 87 IF(ln_wd) THEN 89 ALLOCATE( wd uflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj),STAT=ierr )88 ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj), STAT=ierr ) 90 89 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 91 90 ENDIF 91 ! 92 92 END SUBROUTINE wad_init 93 93 94 94 95 SUBROUTINE wad_lmt( sshb1, sshemp, z2dt ) … … 107 108 ! 108 109 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 109 INTEGER :: zflag ! local scalar110 INTEGER :: jflag ! local scalar 110 111 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars 111 112 REAL(wp) :: zzflxp, zzflxn ! local scalars … … 116 117 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu, zflxv ! local 2D workspace 117 118 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 118 119 119 !!---------------------------------------------------------------------- 120 120 ! … … 131 131 !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 132 132 133 zflag = 0133 jflag = 0 134 134 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 135 135 … … 145 145 ! Horizontal Flux in u and v direction 146 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj m1148 DO ji = 1, jpi m1147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 156 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 157 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 160 161 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells 162 IF(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 158 wdmask(:,:) = 1 159 DO jj = 2, jpj 160 DO ji = 2, jpi 161 162 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE ! we don't care about land cells 163 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 163 164 164 165 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 167 168 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 168 169 169 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 170 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 171 !zdep2 = 0._wp 172 sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 170 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 171 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 172 sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 173 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 174 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 175 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 176 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 177 wdmask(ji,jj) = 0._wp 173 178 END IF 174 179 ENDDO … … 182 187 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 183 188 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 184 185 DO jj = 2, jpjm1 186 DO ji = 2, jpim1 189 jflag = 0 ! flag indicating if any further iterations are needed 190 191 DO jj = 2, jpj 192 DO ji = 2, jpi 187 193 188 wdmask(ji,jj) = 0 189 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 190 IF(bathy(ji,jj) > zdepwd) CYCLE 194 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 195 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 191 196 192 197 ztmp = e1e2t(ji,jj) … … 198 203 199 204 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 200 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) ! this one can be moved out of the loop 201 202 IF(zdep1 > zdep2) THEN 203 zflag = 1 204 wdmask(ji, jj) = 1 205 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 206 207 IF( zdep1 > zdep2 ) THEN 208 wdmask(ji, jj) = 0 205 209 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 206 zcoef = max(zcoef, 0._wp) 210 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 211 ! flag if the limiter has been used but stop flagging if the only 212 ! changes have zeroed the coefficient since further iterations will 213 ! not change anything 214 IF( zcoef > 0._wp ) THEN 215 jflag = 1 216 ELSE 217 zcoef = 0._wp 218 ENDIF 207 219 IF(jk1 > nn_wdit) zcoef = 0._wp 208 220 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 209 221 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 210 222 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 211 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef223 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 212 224 END IF 213 225 END DO ! ji loop … … 217 229 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 218 230 219 IF(lk_mpp) CALL mpp_max(zflag) !max over the global domain 220 221 IF(zflag == 0) EXIT 222 223 zflag = 0 ! flag indicating if any further iteration is needed? 231 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 232 233 IF(jflag == 0) EXIT 234 224 235 END DO ! jk1 loop 225 236 … … 231 242 CALL lbc_lnk( un, 'U', -1. ) 232 243 CALL lbc_lnk( vn, 'V', -1. ) 233 234 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 244 ! 245 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 246 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 247 CALL lbc_lnk( un_b, 'U', -1. ) 248 CALL lbc_lnk( vn_b, 'V', -1. ) 249 250 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 235 251 236 252 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) … … 240 256 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 241 257 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 242 !243 END 244 258 ! 259 ENDIF 260 ! 245 261 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 262 ! 246 263 END SUBROUTINE wad_lmt 264 247 265 248 266 SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt ) … … 260 278 ! 261 279 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 262 INTEGER :: zflag! local scalar280 INTEGER :: jflag ! local scalar 263 281 REAL(wp) :: z2dt 264 282 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars … … 269 287 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! local 2D workspace 270 288 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 271 272 !!---------------------------------------------------------------------- 273 ! 274 289 !!---------------------------------------------------------------------- 290 ! 275 291 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 276 292 … … 284 300 !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 285 301 286 zflag = 0302 jflag = 0 287 303 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 288 304 … … 291 307 zflxp(:,:) = 0._wp 292 308 zflxn(:,:) = 0._wp 293 !zflxu(:,:) = 0._wp294 !zflxv(:,:) = 0._wp295 309 296 310 zwdlmtu(:,:) = 1._wp … … 299 313 ! Horizontal Flux in u and v direction 300 314 301 !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 302 !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 303 304 DO jj = 2, jpjm1 305 DO ji = 2, jpim1 306 307 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells 308 IF(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 315 DO jj = 2, jpj 316 DO ji = 2, jpi 317 318 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells 319 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 309 320 310 321 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 313 324 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 314 325 315 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 316 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 317 !zdep2 = 0._wp 318 sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj) 326 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 327 IF(zdep2 .le. 0._wp) THEN !add more safety, but not necessary 328 sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 329 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 330 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 331 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 332 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 319 333 END IF 320 334 ENDDO … … 328 342 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 329 343 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 330 331 DO jj = 2, jpjm1 332 DO ji = 2, jpim1 344 jflag = 0 ! flag indicating if any further iterations are needed 345 346 DO jj = 2, jpj 347 DO ji = 2, jpi 333 348 334 wdmask(ji,jj) = 0 335 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 336 IF(bathy(ji,jj) > zdepwd) CYCLE 349 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 350 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 337 351 338 352 ztmp = e1e2t(ji,jj) … … 344 358 345 359 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 346 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 ! this one can be moved out of the loop 347 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj) 360 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 348 361 349 362 IF(zdep1 > zdep2) THEN 350 zflag = 1351 !wdmask(ji, jj) = 1352 363 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 353 zcoef = max(zcoef, 0._wp) 364 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 365 ! flag if the limiter has been used but stop flagging if the only 366 ! changes have zeroed the coefficient since further iterations will 367 ! not change anything 368 IF( zcoef > 0._wp ) THEN 369 jflag = 1 370 ELSE 371 zcoef = 0._wp 372 ENDIF 354 373 IF(jk1 > nn_wdit) zcoef = 0._wp 355 374 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 356 375 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 357 376 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 358 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef377 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 359 378 END IF 360 379 END DO ! ji loop … … 364 383 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 365 384 366 IF(lk_mpp) CALL mpp_max(zflag) !max over the global domain 367 368 IF(zflag == 0) EXIT 369 370 zflag = 0 ! flag indicating if any further iteration is needed? 385 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 386 387 IF(jflag == 0) EXIT 388 371 389 END DO ! jk1 loop 372 390 … … 377 395 CALL lbc_lnk( zflxv, 'V', -1. ) 378 396 379 IF( zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'397 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 380 398 381 399 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) … … 385 403 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 386 404 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 387 !405 ! 388 406 END IF 389 407 390 408 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 391 409 END SUBROUTINE wad_lmt_bt 410 411 !!============================================================================== 392 412 END MODULE wet_dry
Note: See TracChangeset
for help on using the changeset viewer.
