Changeset 7421 for branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2016-12-01T18:10:41+01:00 (8 years ago)
- Location:
- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 8 edited
-
dynadv.F90 (modified) (2 diffs)
-
dynhpg.F90 (modified) (18 diffs)
-
dynldf.F90 (modified) (2 diffs)
-
dynspg.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (23 diffs)
-
dynvor.F90 (modified) (13 diffs)
-
dynzdf.F90 (modified) (1 diff)
-
wet_dry.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r6140 r7421 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 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r7412 r7421 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 434 LOGICAL :: ll_tmp1, ll_tmp2, ll_tmp3 ! 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) ) > & 457 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 461 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 462 463 IF(ll_tmp1) THEN 456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) + & 459 & rn_wdmin1 + rn_wdmin2 460 461 IF(ll_tmp1.AND.ll_tmp2) THEN 464 462 zcpx(ji,jj) = 1.0_wp 465 ELSE IF(ll_tmp2) THEN 466 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 467 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 468 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 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) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) / & 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 wduflt(ji,jj) = 1.0_wp 469 469 ELSE 470 470 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp 471 472 END IF 472 473 473 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 474 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 475 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 476 & > rn_wdmin1 + rn_wdmin2 477 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 478 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 479 480 IF(ll_tmp1) THEN 474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) + & 477 & rn_wdmin1 + rn_wdmin2 478 479 IF(ll_tmp1.AND.ll_tmp2) THEN 481 480 zcpy(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp2) 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 )) ) 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) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) / & 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 wdvflt(ji,jj) = 1.0_wp 486 487 ELSE 487 488 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp 488 490 END IF 489 491 END DO 490 492 END DO 491 493 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 492 END IF 494 ENDIF 495 493 496 494 497 ! Surface value … … 507 510 508 511 509 IF( ln_wd) THEN512 IF(ln_wd) THEN 510 513 511 514 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 538 541 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 539 542 540 IF( ln_wd) THEN543 IF(ln_wd) THEN 541 544 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 542 545 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 553 556 ! 554 557 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 555 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )558 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 556 559 ! 557 560 END SUBROUTINE hpg_sco … … 698 701 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 699 702 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 700 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )701 ! 702 ! 703 IF( ln_wd) THEN703 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 704 ! 705 ! 706 IF(ln_wd) THEN 704 707 DO jj = 2, jpjm1 705 708 DO ji = 2, jpim1 706 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 707 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 708 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 709 & > rn_wdmin1 + rn_wdmin2 710 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 711 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) +& 713 & rn_wdmin1 + rn_wdmin2 712 714 713 715 IF(ll_tmp1) THEN 714 716 zcpx(ji,jj) = 1.0_wp 715 717 ELSE IF(ll_tmp2) THEN 716 ! no worries about sshn(ji+1,jj) - sshn(ji,jj) = 0, it won't happen ! here717 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj))&718 & / (sshn(ji+1,jj) - sshn(ji ,jj)))718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 719 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) /& 720 & (sshn(ji+1,jj) - sshn(ji,jj))) 719 721 ELSE 720 722 zcpx(ji,jj) = 0._wp 721 723 END IF 722 724 723 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 724 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 725 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 726 & > rn_wdmin1 + rn_wdmin2 727 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 728 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) & 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) +& 729 & rn_wdmin1 + rn_wdmin2 729 730 730 731 IF(ll_tmp1) THEN 731 732 zcpy(ji,jj) = 1.0_wp 732 733 ELSE IF(ll_tmp2) THEN 733 ! no worries about sshn(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here734 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj))&735 & / (sshn(ji,jj+1) - sshn(ji,jj )))734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 735 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) /& 736 & (sshn(ji,jj+1) - sshn(ji,jj))) 736 737 ELSE 737 738 zcpy(ji,jj) = 0._wp … … 740 741 END DO 741 742 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 742 END IF 743 ENDIF 744 743 745 744 746 IF( kt == nit000 ) THEN … … 911 913 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 912 914 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 913 IF( ln_wd) THEN915 IF(ln_wd) THEN 914 916 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 915 917 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) … … 934 936 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 935 937 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 936 IF( ln_wd) THEN938 IF(ln_wd) THEN 937 939 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 938 940 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 948 950 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 949 951 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 950 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )952 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 951 953 ! 952 954 END SUBROUTINE hpg_djc … … 985 987 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 986 988 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 987 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )989 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 988 990 ! 989 991 IF( kt == nit000 ) THEN … … 998 1000 IF( ln_linssh ) znad = 0._wp 999 1001 1000 IF( ln_wd) THEN1002 IF(ln_wd) THEN 1001 1003 DO jj = 2, jpjm1 1002 1004 DO ji = 2, jpim1 1003 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1004 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 1005 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 1006 & > rn_wdmin1 + rn_wdmin2 1007 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) & 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) +& 1009 & rn_wdmin1 + rn_wdmin2 1009 1010 1010 1011 IF(ll_tmp1) THEN 1011 1012 zcpx(ji,jj) = 1.0_wp 1012 1013 ELSE IF(ll_tmp2) THEN 1013 ! no worries about sshn(ji+1,jj) - sshn(ji,jj) = 0, it won't happen ! here1014 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj))&1015 & / (sshn(ji+1,jj) - sshn(ji ,jj)))1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1015 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) /& 1016 & (sshn(ji+1,jj) - sshn(ji,jj))) 1016 1017 ELSE 1017 1018 zcpx(ji,jj) = 0._wp 1018 1019 END IF 1019 1020 1020 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1021 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 1022 & 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) ) > & 1025 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1) THEN 1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) +& 1025 & rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1.OR.ll_tmp2) THEN 1028 1028 zcpy(ji,jj) = 1.0_wp 1029 1029 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 )))1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 1031 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) /& 1032 & (sshn(ji,jj+1) - sshn(ji,jj))) 1033 1033 ELSE 1034 1034 zcpy(ji,jj) = 0._wp … … 1037 1037 END DO 1038 1038 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 END 1039 ENDIF 1040 1040 1041 1041 ! Clean 3-D work arrays … … 1046 1046 DO jj = 1, jpj 1047 1047 DO ji = 1, jpi 1048 jk = mb athy(ji,jj)1048 jk = mbkt(ji,jj)+1 1049 1049 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 1050 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) … … 1221 1221 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1222 1222 ENDIF 1223 IF( ln_wd) THEN1223 IF(ln_wd) THEN 1224 1224 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 1225 zdpdx2 = zdpdx2 * zcpx(ji,jj) … … 1280 1280 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1281 1281 ENDIF 1282 IF( ln_wd) THEN1282 IF(ln_wd) THEN 1283 1283 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 1284 zdpdy2 = zdpdy2 * zcpy(ji,jj) … … 1295 1295 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1296 1296 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )1297 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1298 1298 ! 1299 1299 END SUBROUTINE hpg_prj -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r6140 r7421 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 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r7412 r7421 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 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7412 r7421 33 33 USE dynvor ! vorticity term 34 34 USE wet_dry ! wetting/drying flux limter 35 USE bdy_ oce , ONLY: ln_bdy35 USE bdy_par ! for lk_bdy 36 36 USE bdytides ! open boundary condition data 37 37 USE bdydyn2d ! open boundary conditions on barotropic variables … … 69 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 70 70 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff /h at F points71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff_f/h at F points 72 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 73 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) … … 156 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 158 159 !!---------------------------------------------------------------------- 159 160 ! … … 167 168 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 168 169 CALL wrk_alloc( jpi,jpj, zhf ) 169 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 170 171 ! 171 172 zmdi=1.e+20 ! missing data indicator for masking … … 219 220 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 220 221 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 221 SELECT CASE( nn_een_e3f ) !* ff /e3 at F-point222 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 222 223 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 223 224 DO jj = 1, jpjm1 … … 225 226 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 226 227 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 227 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff (ji,jj) / zwz(ji,jj)228 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 228 229 END DO 229 230 END DO … … 235 236 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 236 237 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 237 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff (ji,jj) / zwz(ji,jj)238 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 238 239 END DO 239 240 END DO … … 254 255 zwz(:,:) = 0._wp 255 256 zhf(:,:) = 0._wp 256 IF ( .not. ln_sco ) THEN 257 258 !!gm agree the JC comment : this should be done in a much clear way 259 260 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 261 ! Set it to zero for the time being 262 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 263 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 264 ! ENDIF 265 ! zhf(:,:) = gdepw_0(:,:,jk+1) 266 ELSE 267 zhf(:,:) = hbatf(:,:) 268 END IF 269 270 DO jj = 1, jpjm1 271 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 272 END DO 257 258 !!gm assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed 259 !!gm A priori a better value should be something like : 260 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 261 !!gm divided by the sum of the corresponding mask 262 !!gm 263 !! 264 !! IF ( .not. ln_sco ) THEN 265 !! 266 !! !!gm agree the JC comment : this should be done in a much clear way 267 !! 268 !! ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 269 !! ! Set it to zero for the time being 270 !! ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 271 !! ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 272 !! ! ENDIF 273 !! ! zhf(:,:) = gdepw_0(:,:,jk+1) 274 !! ELSE 275 !! zhf(:,:) = hbatf(:,:) 276 !! END IF 277 !! 278 !! DO jj = 1, jpjm1 279 !! zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 280 !! END DO 281 !!gm end 273 282 274 283 DO jk = 1, jpkm1 … … 284 293 END DO 285 294 END DO 286 zwz(:,:) = ff (:,:) * zwz(:,:)295 zwz(:,:) = ff_f(:,:) * zwz(:,:) 287 296 ENDIF 288 297 ENDIF … … 373 382 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 374 383 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 375 DO jj = 2, jpjm1376 DO ji = 2, jpim1377 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &378 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. &379 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) )&380 & > rn_wdmin1 + rn_wdmin2381 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &382 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2383 384 wduflt1(:,:) = 1.0_wp 385 wdvflt1(:,:) = 1.0_wp 386 DO jj = 2, jpjm1 387 DO ji = 2, jpim1 388 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 389 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) & 390 & > rn_wdmin1 + rn_wdmin2 391 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 392 & + rn_wdmin1 + rn_wdmin2 384 393 IF(ll_tmp1) THEN 385 zcpx(ji,jj) = 1.0_wp386 ELSE 387 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen !here388 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) &389 & / (sshn(ji+1,jj) - sshn(ji ,jj)))394 zcpx(ji,jj) = 1.0_wp 395 ELSEIF(ll_tmp2) THEN 396 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 397 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 398 & /(sshn(ji+1,jj) - sshn(ji,jj))) 390 399 ELSE 391 zcpx(ji,jj) = 0._wp 400 zcpx(ji,jj) = 0._wp 401 wduflt1(ji,jj) = 0.0_wp 392 402 END IF 393 394 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 395 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 396 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 397 & > rn_wdmin1 + rn_wdmin2 398 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 399 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 400 403 404 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 405 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) & 406 & > rn_wdmin1 + rn_wdmin2 407 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 408 & + rn_wdmin1 + rn_wdmin2 401 409 IF(ll_tmp1) THEN 402 zcpy(ji,jj)= 1.0_wp403 ELSE 404 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen !here405 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 )))410 zcpy(ji,jj) = 1.0_wp 411 ELSEIF(ll_tmp2) THEN 412 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 413 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 414 & /(sshn(ji,jj+1) - sshn(ji,jj))) 407 415 ELSE 408 zcpy(ji,jj) = 0._wp 409 END IF 410 END DO 416 zcpy(ji,jj) = 0._wp 417 wdvflt1(ji,jj) = 0.0_wp 418 ENDIF 419 420 END DO 411 421 END DO 412 422 423 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 424 413 425 DO jj = 2, jpjm1 414 426 DO ji = 2, jpim1 415 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &416 & * r1_e1u(ji,jj) * zcpx(ji,jj)417 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &418 & * r1_e2v(ji,jj) * zcpy(ji,jj)427 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 428 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 429 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 430 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 419 431 END DO 420 432 END DO … … 563 575 ENDIF 564 576 577 IF( ln_wd ) THEN !preserve the positivity of water depth 578 !ssh[b,n,a] should have already been processed for this 579 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - ht_0(:,:)) 580 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - ht_0(:,:)) 581 ENDIF 565 582 ! 566 583 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 598 615 ! ! ------------------ 599 616 ! Update only tidal forcing at open boundaries 600 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 601 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 617 #if defined key_tide 618 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 619 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 620 #endif 602 621 ! 603 622 ! Set extrapolation coefficients for predictor step: … … 635 654 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 636 655 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 656 IF( ln_wd ) THEN 657 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 658 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 659 END IF 637 660 ELSE 638 661 zhup2_e (:,:) = hu_n(:,:) … … 686 709 END DO 687 710 END DO 688 689 711 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 690 712 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - ht_0(:,:)) 691 713 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 692 714 715 #if defined key_bdy 693 716 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 694 IF( l n_bdy ) CALL bdy_ssh( ssha_e )695 717 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 718 #endif 696 719 #if defined key_agrif 697 720 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) … … 734 757 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 735 758 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 736 737 759 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 760 wduflt1(:,:) = 1._wp 761 wdvflt1(:,:) = 1._wp 738 762 DO jj = 2, jpjm1 739 DO ji = 2, jpim1 740 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 741 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 742 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 743 & > rn_wdmin1 + rn_wdmin2 744 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 745 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 746 747 IF(ll_tmp1) THEN 748 zcpx(ji,jj) = 1.0_wp 749 ELSE IF(ll_tmp2) THEN 750 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 751 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 752 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 753 ELSE 754 zcpx(ji,jj) = 0._wp 755 END IF 756 757 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 758 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 759 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 760 & > rn_wdmin1 + rn_wdmin2 761 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 762 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 763 764 IF(ll_tmp1) THEN 765 zcpy(ji,jj) = 1.0_wp 766 ELSE IF(ll_tmp2) THEN 767 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 768 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 769 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 770 ELSE 771 zcpy(ji,jj) = 0._wp 772 END IF 763 DO ji = 2, jpim1 764 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -ht_0(ji,jj), -ht_0(ji+1,jj) ) & 765 & .AND. MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 766 & > rn_wdmin1 + rn_wdmin2 767 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -ht_0(ji,jj), -ht_0(ji+1,jj) ) & 768 & + rn_wdmin1 + rn_wdmin2 769 IF(ll_tmp1) THEN 770 zcpx(ji,jj) = 1._wp 771 ELSE IF(ll_tmp2) THEN 772 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 773 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 774 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 775 ELSE 776 zcpx(ji,jj) = 0._wp 777 wduflt1(ji,jj) = 0._wp 778 END IF 779 780 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -ht_0(ji,jj), -ht_0(ji,jj+1) ) & 781 & .AND. MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 782 & > rn_wdmin1 + rn_wdmin2 783 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -ht_0(ji,jj), -ht_0(ji,jj+1) ) & 784 & + rn_wdmin1 + rn_wdmin2 785 IF(ll_tmp1) THEN 786 zcpy(ji,jj) = 1._wp 787 ELSE IF(ll_tmp2) THEN 788 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 789 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 790 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 791 ELSE 792 zcpy(ji,jj) = 0._wp 793 wdvflt1(ji,jj) = 0._wp 794 END IF 773 795 END DO 774 END DO 775 END IF 796 END DO 797 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 798 ENDIF 776 799 ! 777 800 ! Compute associated depths at U and V points: … … 791 814 END DO 792 815 816 IF( ln_wd ) THEN 817 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 818 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 819 END IF 820 793 821 ENDIF 794 822 ! … … 841 869 ! 842 870 ! Add tidal astronomical forcing if defined 843 IF ( l n_tide.AND.ln_tide_pot ) THEN871 IF ( lk_tide.AND.ln_tide_pot ) THEN 844 872 DO jj = 2, jpjm1 845 873 DO ji = fs_2, fs_jpim1 ! vector opt. … … 868 896 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 869 897 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 870 zwx(ji,jj) = zu_spg * zcpx(ji,jj) * wdmask(ji,jj) * wdmask(ji+1, jj)871 zwy(ji,jj) = zv_spg * zcpy(ji,jj) * wdmask(ji,jj) * wdmask(ji, jj+1)898 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 899 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 872 900 END DO 873 901 END DO … … 907 935 DO ji = fs_2, fs_jpim1 ! vector opt. 908 936 909 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 910 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 937 IF( ln_wd ) THEN 938 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 939 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 940 ELSE 941 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 942 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 943 END IF 911 944 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 912 945 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 928 961 ! 929 962 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 930 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 931 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 963 IF( ln_wd ) THEN 964 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 965 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 966 ELSE 967 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 968 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 969 END IF 932 970 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 933 971 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 937 975 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 938 976 ! 977 #if defined key_bdy 939 978 ! ! open boundaries 940 IF( l n_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )941 979 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 980 #endif 942 981 #if defined key_agrif 943 982 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif … … 993 1032 ! 994 1033 ! Update barotropic trend: 995 IF(ln_wd) THEN 996 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 997 DO jk=1,jpkm1 998 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b * wdmask(:,:) 999 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1000 END DO 1001 ELSE 1002 ! At this stage, ssha has been corrected: compute new depths at velocity points 1003 DO jj = 1, jpjm1 1004 DO ji = 1, jpim1 ! NO Vector Opt. 1005 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1006 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1007 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1008 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1009 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1010 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1011 END DO 1012 END DO 1013 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1014 ! 1015 DO jk=1,jpkm1 1016 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1017 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1018 END DO 1019 ! Save barotropic velocities not transport: 1020 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1021 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1022 ENDIF 1034 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1035 DO jk=1,jpkm1 1036 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1037 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1038 END DO 1023 1039 ELSE 1024 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1025 DO jk=1,jpkm1 1026 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1027 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1028 END DO 1029 ELSE 1030 ! At this stage, ssha has been corrected: compute new depths at velocity points 1031 DO jj = 1, jpjm1 1032 DO ji = 1, jpim1 ! NO Vector Opt. 1033 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1034 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1035 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1036 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1037 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1038 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1039 END DO 1040 END DO 1041 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1042 ! 1043 DO jk=1,jpkm1 1044 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1045 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1046 END DO 1047 ! Save barotropic velocities not transport: 1048 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1049 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1050 ENDIF 1051 1052 END IF 1040 ! At this stage, ssha has been corrected: compute new depths at velocity points 1041 DO jj = 1, jpjm1 1042 DO ji = 1, jpim1 ! NO Vector Opt. 1043 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1044 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1045 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1046 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1047 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1048 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1049 END DO 1050 END DO 1051 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1052 ! 1053 DO jk=1,jpkm1 1054 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1055 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1056 END DO 1057 ! Save barotropic velocities not transport: 1058 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1059 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1060 ENDIF 1053 1061 ! 1054 1062 DO jk = 1, jpkm1 … … 1086 1094 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1087 1095 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy )1096 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 1089 1097 ! 1090 1098 IF ( ln_diatmb ) THEN -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r6140 r7421 237 237 SELECT CASE( kvor ) !== vorticity considered ==! 238 238 CASE ( np_COR ) !* Coriolis (planetary vorticity) 239 zwz(:,:) = ff (:,:)239 zwz(:,:) = ff_f(:,:) 240 240 CASE ( np_RVO ) !* relative vorticity 241 241 DO jj = 1, jpjm1 … … 256 256 DO jj = 1, jpjm1 257 257 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) ) &260 & * r1_e1e2f(ji,jj)258 zwz(ji,jj) = ff_f(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) ) & 260 & * r1_e1e2f(ji,jj) 261 261 END DO 262 262 END DO … … 264 264 DO jj = 1, jpjm1 265 265 DO ji = 1, fs_jpim1 ! vector opt. 266 zwz(ji,jj) = ff (ji,jj)&266 zwz(ji,jj) = ff_f(ji,jj) & 267 267 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 268 268 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & … … 357 357 SELECT CASE( kvor ) !== vorticity considered ==! 358 358 CASE ( np_COR ) !* Coriolis (planetary vorticity) 359 zwz(:,:) = ff (:,:)359 zwz(:,:) = ff_f(:,:) 360 360 CASE ( np_RVO ) !* relative vorticity 361 361 DO jj = 1, jpjm1 … … 376 376 DO jj = 1, jpjm1 377 377 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) ) &380 & * r1_e1e2f(ji,jj)378 zwz(ji,jj) = ff_f(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) ) & 380 & * r1_e1e2f(ji,jj) 381 381 END DO 382 382 END DO … … 384 384 DO jj = 1, jpjm1 385 385 DO ji = 1, fs_jpim1 ! vector opt. 386 zwz(ji,jj) = ff (ji,jj)&386 zwz(ji,jj) = ff_f(ji,jj) & 387 387 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 388 388 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & … … 506 506 DO jj = 1, jpjm1 507 507 DO ji = 1, fs_jpim1 ! vector opt. 508 zwz(ji,jj) = ff (ji,jj) * z1_e3f(ji,jj)508 zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj) 509 509 END DO 510 510 END DO … … 528 528 DO jj = 1, jpjm1 529 529 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) ) &532 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)530 zwz(ji,jj) = ( ff_f(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) ) & 532 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 533 533 END DO 534 534 END DO … … 536 536 DO jj = 1, jpjm1 537 537 DO ji = 1, fs_jpim1 ! vector opt. 538 zwz(ji,jj) = ( ff (ji,jj)&538 zwz(ji,jj) = ( ff_f(ji,jj) & 539 539 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 540 540 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & … … 626 626 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 627 627 WRITE(numout,*) '~~~~~~~~~~~~' 628 WRITE(numout,*) ' 629 WRITE(numout,*) ' 630 WRITE(numout,*) ' 631 WRITE(numout,*) ' 632 WRITE(numout,*) ' 633 WRITE(numout,*) ' 634 WRITE(numout,*) ' 628 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 629 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 630 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 631 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 632 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 633 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 634 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 635 635 ENDIF 636 636 … … 639 639 ! at angles with three ocean points and one land point 640 640 IF(lwp) WRITE(numout,*) 641 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T)ln_vorlat = ', ln_vorlat641 IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 642 642 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 643 643 DO jk = 1, jpk … … 666 666 ncor = np_COR 667 667 IF( ln_dynadv_vec ) THEN 668 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'668 IF(lwp) WRITE(numout,*) ' ===>> Vector form advection : vorticity = Coriolis + relative vorticity' 669 669 nrvm = np_RVO ! relative vorticity 670 670 ntot = np_CRV ! relative + planetary vorticity 671 671 ELSE 672 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'672 IF(lwp) WRITE(numout,*) ' ===>> Flux form advection : vorticity = Coriolis + metric term' 673 673 nrvm = np_MET ! metric term 674 674 ntot = np_CME ! Coriolis + metric term … … 677 677 IF(lwp) THEN ! Print the choice 678 678 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'679 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' ===>> energy conserving scheme' 680 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' ===>> enstrophy conserving scheme' 681 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' ===>> mixed enstrophy/energy conserving scheme' 682 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' ===>> energy and enstrophy conserving scheme' 683 683 ENDIF 684 684 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r6140 r7421 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 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r7412 r7421 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 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter 35 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 36 35 … … 45 44 PUBLIC wad_lmt ! routine called by sshwzv.F90 46 45 PUBLIC wad_lmt_bt ! routine called by dynspg_ts.F90 47 PUBLIC wad_istate ! routine called by istate.F90 and domvvl.F9048 46 49 47 !! * Substitutions … … 77 75 WRITE(numout,*) 78 76 WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read' 79 WRITE(numout,*) '~~~~~~~ 77 WRITE(numout,*) '~~~~~~~~' 80 78 WRITE(numout,*) ' Namelist namwad' 81 79 WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd … … 84 82 WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld 85 83 WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit 86 87 84 ENDIF 85 ! 88 86 IF(ln_wd) THEN 89 ALLOCATE( wd mask(jpi,jpj), STAT=ierr )87 ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr ) 90 88 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 91 89 ENDIF 90 ! 92 91 END SUBROUTINE wad_init 92 93 93 94 94 SUBROUTINE wad_lmt( sshb1, sshemp, z2dt ) … … 116 116 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu, zflxv ! local 2D workspace 117 117 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 118 119 118 !!---------------------------------------------------------------------- 120 119 ! … … 124 123 IF(ln_wd) THEN 125 124 126 CALL wrk_alloc( jpi, jpj,zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )127 CALL wrk_alloc( jpi, jpj,zwdlmtu, zwdlmtv)125 CALL wrk_alloc( jpi,jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 126 CALL wrk_alloc( jpi,jpj, zwdlmtu, zwdlmtv) 128 127 ! 129 128 … … 145 144 ! Horizontal Flux in u and v direction 146 145 DO jk = 1, jpkm1 147 DO jj = 1, jpj 148 DO ji = 1, jpi 146 DO jj = 1, jpjm1 147 DO ji = 1, jpim1 149 148 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 149 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 156 155 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 156 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(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 157 DO jj = 2, jpjm1 158 DO ji = 2, jpim1 159 160 IF( tmask(ji,jj,1) == 0._wp ) CYCLE ! we don't care about land cells 161 IF( ht_0 (ji,jj) > zdepwd ) CYCLE ! and cells which will unlikely go dried out 164 162 165 163 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 168 166 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 169 167 170 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1171 IF(zdep2 .le.0._wp) THEN !add more safty, but not necessary168 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 169 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 172 170 !zdep2 = 0._wp 173 sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 174 wdmask(ji,jj) = 0._wp 171 sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 175 172 END IF 176 173 ENDDO … … 185 182 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 186 183 187 DO jj = 2, jpj 188 DO ji = 2, jpi 184 DO jj = 2, jpjm1 185 DO ji = 2, jpim1 189 186 190 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 191 IF(bathy(ji,jj) > zdepwd) CYCLE 187 wdmask(ji,jj) = 0 188 IF( tmask(ji,jj,1) < 0.5_wp) CYCLE 189 IF( ht_0(ji,jj) > zdepwd) CYCLE 192 190 193 191 ztmp = e1e2t(ji,jj) … … 199 197 200 198 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 201 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) ! this one can be moved out of the loop199 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) ! this one can be moved out of the loop 202 200 203 201 IF(zdep1 > zdep2) THEN 204 202 zflag = 1 205 wdmask(ji, jj) = 0203 wdmask(ji, jj) = 1 206 204 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 207 205 zcoef = max(zcoef, 0._wp) … … 210 208 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 211 209 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 212 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji ,jj-1) = zcoef210 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 213 211 END IF 214 212 END DO ! ji loop … … 232 230 CALL lbc_lnk( un, 'U', -1. ) 233 231 CALL lbc_lnk( vn, 'V', -1. ) 234 !235 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)236 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)237 CALL lbc_lnk( un_b, 'U', -1. )238 CALL lbc_lnk( vn_b, 'V', -1. )239 232 240 233 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' … … 246 239 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 247 240 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 248 !249 END 250 241 ! 242 ENDIF 243 ! 251 244 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 245 ! 252 246 END SUBROUTINE wad_lmt 247 253 248 254 249 SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt ) … … 275 270 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! local 2D workspace 276 271 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 277 278 !!---------------------------------------------------------------------- 279 ! 280 272 !!---------------------------------------------------------------------- 273 ! 281 274 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 282 275 … … 297 290 zflxp(:,:) = 0._wp 298 291 zflxn(:,:) = 0._wp 292 !zflxu(:,:) = 0._wp 293 !zflxv(:,:) = 0._wp 299 294 300 295 zwdlmtu(:,:) = 1._wp … … 303 298 ! Horizontal Flux in u and v direction 304 299 305 DO jj = 2, jpj 306 DO ji = 2, jpi 307 308 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells 309 IF(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 300 !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 301 !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 302 303 DO jj = 2, jpjm1 304 DO ji = 2, jpim1 305 306 IF(tmask(ji,jj,1) < 0.5_wp) CYCLE ! we don't care about land cells 307 IF(ht_0 (ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 310 308 311 309 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 314 312 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 315 313 316 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 314 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 315 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 316 !zdep2 = 0._wp 317 sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 318 END IF 317 319 ENDDO 318 320 END DO … … 326 328 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 327 329 328 DO jj = 2, jpj 329 DO ji = 2, jpi 330 DO jj = 2, jpjm1 331 DO ji = 2, jpim1 330 332 331 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 332 IF(bathy(ji,jj) > zdepwd) CYCLE 333 wdmask(ji,jj) = 0 334 IF(tmask(ji,jj,1) < 0.5_wp) CYCLE 335 IF(ht_0 (ji,jj) > zdepwd) CYCLE 333 336 334 337 ztmp = e1e2t(ji,jj) … … 340 343 341 344 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 342 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 ! this one can be moved out of the loop345 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 ! this one can be moved out of the loop 343 346 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj) 344 347 345 348 IF(zdep1 > zdep2) THEN 346 349 zflag = 1 350 !wdmask(ji, jj) = 1 347 351 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 348 352 zcoef = max(zcoef, 0._wp) … … 351 355 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 352 356 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 353 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji ,jj-1) = zcoef357 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 354 358 END IF 355 359 END DO ! ji loop … … 374 378 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 375 379 380 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 381 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 376 382 ! 377 383 ! 378 384 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 379 385 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 380 !386 ! 381 387 END IF 382 388 ! 383 389 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 390 ! 384 391 END SUBROUTINE wad_lmt_bt 385 392 386 SUBROUTINE wad_istate 387 !!---------------------------------------------------------------------- 388 !! *** ROUTINE wad_istate *** 389 !! 390 !! ** Purpose : Initialization of the dynamics and tracers for WAD test 391 !! configurations (channels or bowls with initial ssh gradients) 392 !! 393 !! ** Method : - set temperature field 394 !! - set salinity field 395 !! - set ssh slope (needs to be repeated in domvvl_rst_init to 396 !! set vertical metrics ) 397 !!---------------------------------------------------------------------- 398 ! 399 INTEGER :: ji, jj ! dummy loop indices 400 REAL(wp) :: zi, zj 401 !!---------------------------------------------------------------------- 402 ! 403 ! Uniform T & S in all test cases 404 tsn(:,:,:,jp_tem) = 10._wp 405 tsb(:,:,:,jp_tem) = 10._wp 406 tsn(:,:,:,jp_sal) = 35._wp 407 tsb(:,:,:,jp_sal) = 35._wp 408 SELECT CASE ( jp_cfg ) 409 ! ! ==================== 410 CASE ( 1 ) ! WAD 1 configuration 411 ! ! ==================== 412 ! 413 IF(lwp) WRITE(numout,*) 414 IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope' 415 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 416 ! 417 do ji = 1,jpi 418 sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 419 end do 420 ! ! ==================== 421 CASE ( 2 ) ! WAD 2 configuration 422 ! ! ==================== 423 ! 424 IF(lwp) WRITE(numout,*) 425 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope' 426 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 427 ! 428 do ji = 1,jpi 429 sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 430 end do 431 ! ! ==================== 432 CASE ( 3 ) ! WAD 3 configuration 433 ! ! ==================== 434 ! 435 IF(lwp) WRITE(numout,*) 436 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 437 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 438 ! 439 do ji = 1,jpi 440 sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 441 end do 442 443 ! 444 ! ! ==================== 445 CASE ( 4 ) ! WAD 4 configuration 446 ! ! ==================== 447 ! 448 IF(lwp) WRITE(numout,*) 449 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 450 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 451 ! 452 DO ji = 1, jpi 453 zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 ) 454 DO jj = 1, jpj 455 zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 ) 456 sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj 457 END DO 458 END DO 459 460 ! 461 ! ! =========================== 462 CASE ( 5 ) ! WAD 5 configuration 463 ! ! ==================== 464 ! 465 IF(lwp) WRITE(numout,*) 466 IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf' 467 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 468 ! 469 ! Needed rn_wdmin2 increased to 0.01 for this case? 470 do ji = 1,jpi 471 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 472 end do 473 474 ! 475 ! ! =========================== 476 CASE ( 6 ) ! WAD 6 configuration 477 ! ! ==================== 478 ! 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 482 ! 483 do ji = 1,jpi 484 !6a 485 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 486 !Some variations in initial slope that have been tested 487 !6b 488 !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 489 !6c 490 !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 491 !6d 492 !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 493 end do 494 495 ! 496 ! ! =========================== 497 CASE DEFAULT ! NONE existing configuration 498 ! ! =========================== 499 WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 500 ! 501 CALL ctl_stop( ctmp1 ) 502 ! 503 END SELECT 504 ! 505 ! Apply minimum wetdepth criterion 506 ! 507 do jj = 1,jpj 508 do ji = 1,jpi 509 IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN 510 sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) ) 511 ENDIF 512 end do 513 end do 514 sshb = sshn 515 ssha = sshn 516 ! 517 END SUBROUTINE wad_istate 518 519 !!===================================================================== 393 !!============================================================================== 520 394 END MODULE wet_dry
Note: See TracChangeset
for help on using the changeset viewer.
