- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5506 r6808 2 2 !!============================================================================== 3 3 !! *** MODULE domzgr *** 4 !! Ocean initialization : domain initialization4 !! Ocean domain : definition of the vertical coordinate system 5 5 !!============================================================================== 6 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate … … 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 20 21 !!---------------------------------------------------------------------- 21 22 … … 36 37 USE oce ! ocean variables 37 38 USE dom_oce ! ocean domain 39 USE wet_dry ! wetting and drying 38 40 USE closea ! closed seas 39 41 USE c1d ! 1D vertical configuration 42 ! 40 43 USE in_out_manager ! I/O manager 41 44 USE iom ! I/O library … … 73 76 74 77 !! * Substitutions 75 # include "domzgr_substitute.h90"76 78 # include "vectopt_loop_substitute.h90" 77 79 !!---------------------------------------------------------------------- … … 102 104 INTEGER :: ios 103 105 ! 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 105 107 !!---------------------------------------------------------------------- 106 108 ! … … 120 122 WRITE(numout,*) 'dom_zgr : vertical coordinate' 121 123 WRITE(numout,*) '~~~~~~~' 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 125 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 129 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 127 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 133 129 134 ioptio = 0 ! Check Vertical coordinate options … … 155 160 ! 156 161 IF( nprint == 1 .AND. lwp ) THEN 157 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )162 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 158 163 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 159 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )160 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &161 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &162 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &163 & ' w ', MINVAL(e3w_0(:,:,:) )164 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 165 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 166 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 167 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 168 & ' w ', MINVAL( e3w_0(:,:,:) ) 164 169 165 170 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 166 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )167 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &168 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &169 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&170 & ' w ', MAXVAL(e3w_0(:,:,:) )171 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 172 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 173 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 174 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 175 & ' w ', MAXVAL( e3w_0(:,:,:) ) 171 176 ENDIF 172 177 ! … … 219 224 & ppsur == pp_to_be_computed ) THEN 220 225 ! 226 #if defined key_agrif 227 za1 = ( ppdzmin - pphmax / FLOAT(jpkdta-1) ) & 228 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * ( LOG( COSH( (jpkdta - ppkth) / ppacr) )& 229 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 230 #else 221 231 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 222 232 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 223 233 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 234 #endif 224 235 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 225 236 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) … … 236 247 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 237 248 WRITE(numout,*) ' Total depth :', zhmax 249 #if defined key_agrif 250 WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1) 251 #else 238 252 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 253 #endif 239 254 ELSE 240 255 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN … … 260 275 ! Reference z-coordinate (depth - scale factor at T- and W-points) 261 276 ! ====================== 262 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 277 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 278 #if defined key_agrif 279 za1 = zhmax / FLOAT(jpkdta-1) 280 #else 263 281 za1 = zhmax / FLOAT(jpk-1) 282 #endif 264 283 DO jk = 1, jpk 265 284 zw = FLOAT( jk ) … … 365 384 !! - bathy : meter bathymetry (in meters) 366 385 !!---------------------------------------------------------------------- 367 INTEGER :: ji, jj, j l, jk ! dummy loop indices386 INTEGER :: ji, jj, jk ! dummy loop indices 368 387 INTEGER :: inum ! temporary logical unit 369 388 INTEGER :: ierror ! error flag … … 487 506 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 488 507 ! ! ===================== 489 IF( nn_cla == 0 ) THEN 490 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 491 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 492 DO ji = mi0(ii0), mi1(ii1) 493 DO jj = mj0(ij0), mj1(ij1) 494 mbathy(ji,jj) = 15 495 END DO 508 ! 509 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 510 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 511 DO ji = mi0(ii0), mi1(ii1) 512 DO jj = mj0(ij0), mj1(ij1) 513 mbathy(ji,jj) = 15 496 514 END DO 497 IF(lwp) WRITE(numout,*)498 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0499 !500 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open501 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)502 DO ji = mi0(ii0), mi1(ii1)503 DO jj = mj0(ij0), mj1(ij1)504 mbathy(ji,jj) = 12505 END DO515 END DO 516 IF(lwp) WRITE(numout,*) 517 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 518 ! 519 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 520 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 521 DO ji = mi0(ii0), mi1(ii1) 522 DO jj = mj0(ij0), mj1(ij1) 523 mbathy(ji,jj) = 12 506 524 END DO 507 IF(lwp) WRITE(numout,*)508 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0509 ENDIF525 END DO 526 IF(lwp) WRITE(numout,*) 527 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 510 528 ! 511 529 ENDIF … … 528 546 CALL iom_close( inum ) 529 547 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 548 549 ! set grounded point to 0 550 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 551 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 552 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 553 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 554 END WHERE 530 555 END IF 531 556 ! 532 557 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 533 ! 534 IF( nn_cla == 0 ) THEN 535 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 536 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 537 DO ji = mi0(ii0), mi1(ii1) 538 DO jj = mj0(ij0), mj1(ij1) 539 bathy(ji,jj) = 284._wp 540 END DO 558 ! 559 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 560 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 561 DO ji = mi0(ii0), mi1(ii1) 562 DO jj = mj0(ij0), mj1(ij1) 563 bathy(ji,jj) = 284._wp 541 564 END DO 542 IF(lwp) WRITE(numout,*)543 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0544 !545 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open546 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)547 DO ji = mi0(ii0), mi1(ii1)548 DO jj = mj0(ij0), mj1(ij1)549 bathy(ji,jj) = 137._wp550 END DO565 END DO 566 IF(lwp) WRITE(numout,*) 567 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 568 ! 569 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 570 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 571 DO ji = mi0(ii0), mi1(ii1) 572 DO jj = mj0(ij0), mj1(ij1) 573 bathy(ji,jj) = 137._wp 551 574 END DO 552 IF(lwp) WRITE(numout,*)553 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0554 ENDIF575 END DO 576 IF(lwp) WRITE(numout,*) 577 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 555 578 ! 556 579 ENDIF … … 567 590 ! 568 591 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 569 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin570 IF ( ln_isfcav ) THEN571 WHERE (bathy == risfdep)572 bathy = 0.0_wp ; risfdep = 0.0_wp573 END WHERE574 END IF575 ! end patch576 592 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 577 593 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 663 679 !! - update bathy : meter bathymetry (in meters) 664 680 !!---------------------------------------------------------------------- 665 !!666 681 INTEGER :: ji, jj, jl ! dummy loop indices 667 682 INTEGER :: icompt, ibtest, ikmax ! temporary integers 668 683 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 669 670 684 !!---------------------------------------------------------------------- 671 685 ! … … 764 778 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 765 779 ENDIF 766 767 IF( lwp .AND. nprint == 1 ) THEN ! control print768 WRITE(numout,*)769 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '770 WRITE(numout,*) ' ------------------'771 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )772 WRITE(numout,*)773 ENDIF774 780 ! 775 781 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 792 798 !! (min value = 1 over land) 793 799 !!---------------------------------------------------------------------- 794 !!795 800 INTEGER :: ji, jj ! dummy loop indices 796 801 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 824 829 END SUBROUTINE zgr_bot_level 825 830 826 SUBROUTINE zgr_top_level 827 !!---------------------------------------------------------------------- 828 !! *** ROUTINE zgr_bot_level *** 831 832 SUBROUTINE zgr_top_level 833 !!---------------------------------------------------------------------- 834 !! *** ROUTINE zgr_top_level *** 829 835 !! 830 836 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 836 842 !! (min value = 1) 837 843 !!---------------------------------------------------------------------- 838 !!839 844 INTEGER :: ji, jj ! dummy loop indices 840 845 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 870 875 END SUBROUTINE zgr_top_level 871 876 877 872 878 SUBROUTINE zgr_zco 873 879 !!---------------------------------------------------------------------- 874 880 !! *** ROUTINE zgr_zco *** 875 881 !! 876 !! ** Purpose : define the z-coordinate system882 !! ** Purpose : define the reference z-coordinate system 877 883 !! 878 884 !! ** Method : set 3D coord. arrays to reference 1D array … … 884 890 ! 885 891 DO jk = 1, jpk 886 gdept_0 887 gdepw_0 888 gde p3w_0(:,:,jk) = gdepw_1d(jk)889 e3t_0 890 e3u_0 891 e3v_0 892 e3f_0 893 e3w_0 894 e3uw_0 895 e3vw_0 892 gdept_0(:,:,jk) = gdept_1d(jk) 893 gdepw_0(:,:,jk) = gdepw_1d(jk) 894 gde3w_0(:,:,jk) = gdepw_1d(jk) 895 e3t_0 (:,:,jk) = e3t_1d (jk) 896 e3u_0 (:,:,jk) = e3t_1d (jk) 897 e3v_0 (:,:,jk) = e3t_1d (jk) 898 e3f_0 (:,:,jk) = e3t_1d (jk) 899 e3w_0 (:,:,jk) = e3w_1d (jk) 900 e3uw_0 (:,:,jk) = e3w_1d (jk) 901 e3vw_0 (:,:,jk) = e3w_1d (jk) 896 902 END DO 897 903 ! … … 906 912 !! 907 913 !! ** Purpose : the depth and vertical scale factor in partial step 908 !! z-coordinate case914 !! reference z-coordinate case 909 915 !! 910 916 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 946 952 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 947 953 !!---------------------------------------------------------------------- 948 !!949 954 INTEGER :: ji, jj, jk ! dummy loop indices 950 955 INTEGER :: ik, it, ikb, ikt ! temporary integers 951 LOGICAL :: ll_print ! Allow control print for debugging952 956 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 953 957 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 954 REAL(wp) :: zmax ! Maximum depth955 958 REAL(wp) :: zdiff ! temporary scalar 956 REAL(wp) :: z refdep! temporary scalar959 REAL(wp) :: zmax ! temporary scalar 957 960 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 958 961 !!--------------------------------------------------------------------- … … 960 963 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 961 964 ! 962 CALL wrk_alloc( jpi, jpj, jpk,zprt )965 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 963 966 ! 964 967 IF(lwp) WRITE(numout,*) … … 966 969 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 967 970 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 968 969 ll_print = .FALSE. ! Local variable for debugging970 971 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth972 WRITE(numout,*)973 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'974 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )975 ENDIF976 977 971 978 972 ! bathymetry in level (from bathy_meter) … … 993 987 END DO 994 988 995 IF ( ln_isfcav ) CALL zgr_isf996 997 989 ! Scale factors and depth at T- and W-points 998 990 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1002 994 e3w_0 (:,:,jk) = e3w_1d (jk) 1003 995 END DO 996 997 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 998 IF ( ln_isfcav ) CALL zgr_isf 999 1000 ! Scale factors and depth at T- and W-points 1001 IF ( .NOT. ln_isfcav ) THEN 1002 DO jj = 1, jpj 1003 DO ji = 1, jpi 1004 ik = mbathy(ji,jj) 1005 IF( ik > 0 ) THEN ! ocean point only 1006 ! max ocean level case 1007 IF( ik == jpkm1 ) THEN 1008 zdepwp = bathy(ji,jj) 1009 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1010 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1011 e3t_0(ji,jj,ik ) = ze3tp 1012 e3t_0(ji,jj,ik+1) = ze3tp 1013 e3w_0(ji,jj,ik ) = ze3wp 1014 e3w_0(ji,jj,ik+1) = ze3tp 1015 gdepw_0(ji,jj,ik+1) = zdepwp 1016 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1017 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1018 ! 1019 ELSE ! standard case 1020 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1021 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1022 ENDIF 1023 !gm Bug? check the gdepw_1d 1024 ! ... on ik 1025 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1026 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1027 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1028 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1029 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1030 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1031 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1032 ! ... on ik+1 1033 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1034 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1035 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1036 ENDIF 1037 ENDIF 1038 END DO 1039 END DO 1040 ! 1041 it = 0 1042 DO jj = 1, jpj 1043 DO ji = 1, jpi 1044 ik = mbathy(ji,jj) 1045 IF( ik > 0 ) THEN ! ocean point only 1046 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1047 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1048 ! test 1049 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1050 IF( zdiff <= 0._wp .AND. lwp ) THEN 1051 it = it + 1 1052 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1053 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1054 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1055 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1056 ENDIF 1057 ENDIF 1058 END DO 1059 END DO 1060 END IF 1061 ! 1062 ! Scale factors and depth at U-, V-, UW and VW-points 1063 DO jk = 1, jpk ! initialisation to z-scale factors 1064 e3u_0 (:,:,jk) = e3t_1d(jk) 1065 e3v_0 (:,:,jk) = e3t_1d(jk) 1066 e3uw_0(:,:,jk) = e3w_1d(jk) 1067 e3vw_0(:,:,jk) = e3w_1d(jk) 1068 END DO 1069 1070 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1071 DO jj = 1, jpjm1 1072 DO ji = 1, fs_jpim1 ! vector opt. 1073 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1074 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1075 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1076 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1077 END DO 1078 END DO 1079 END DO 1080 IF ( ln_isfcav ) THEN 1081 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1082 DO jj = 2, jpjm1 1083 DO ji = 2, fs_jpim1 ! vector opt. 1084 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1085 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1086 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1087 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1088 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1089 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1090 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1091 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1092 END DO 1093 END DO 1094 END IF 1095 1096 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1097 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1098 ! 1099 1100 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1101 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1102 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1103 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1104 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1105 END DO 1106 1107 ! Scale factor at F-point 1108 DO jk = 1, jpk ! initialisation to z-scale factors 1109 e3f_0(:,:,jk) = e3t_1d(jk) 1110 END DO 1111 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1112 DO jj = 1, jpjm1 1113 DO ji = 1, fs_jpim1 ! vector opt. 1114 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1115 END DO 1116 END DO 1117 END DO 1118 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1119 ! 1120 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1121 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1122 END DO 1123 !!gm bug ? : must be a do loop with mj0,mj1 1004 1124 ! 1125 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1126 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1127 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1128 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1129 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1130 1131 ! Control of the sign 1132 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1133 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1134 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1135 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1136 1137 ! Compute gde3w_0 (vertical sum of e3w) 1138 IF ( ln_isfcav ) THEN ! if cavity 1139 WHERE( misfdep == 0 ) misfdep = 1 1140 DO jj = 1,jpj 1141 DO ji = 1,jpi 1142 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1143 DO jk = 2, misfdep(ji,jj) 1144 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1145 END DO 1146 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1147 DO jk = misfdep(ji,jj) + 1, jpk 1148 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1149 END DO 1150 END DO 1151 END DO 1152 ELSE ! no cavity 1153 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1154 DO jk = 2, jpk 1155 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1156 END DO 1157 END IF 1158 ! 1159 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1160 ! 1161 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1162 ! 1163 END SUBROUTINE zgr_zps 1164 1165 1166 SUBROUTINE zgr_isf 1167 !!---------------------------------------------------------------------- 1168 !! *** ROUTINE zgr_isf *** 1169 !! 1170 !! ** Purpose : check the bathymetry in levels 1171 !! 1172 !! ** Method : THe water column have to contained at least 2 cells 1173 !! Bathymetry and isfdraft are modified (dig/close) to respect 1174 !! this criterion. 1175 !! 1176 !! ** Action : - test compatibility between isfdraft and bathy 1177 !! - bathy and isfdraft are modified 1178 !!---------------------------------------------------------------------- 1179 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1180 INTEGER :: ik, it ! temporary integers 1181 INTEGER :: icompt, ibtest ! (ISF) 1182 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1183 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1184 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1185 REAL(wp) :: zmax ! Maximum and minimum depth 1186 REAL(wp) :: zbathydiff ! isf temporary scalar 1187 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1188 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1189 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1190 REAL(wp) :: zdiff ! temporary scalar 1191 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1192 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1193 !!--------------------------------------------------------------------- 1194 ! 1195 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1196 ! 1197 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1198 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1199 1200 ! (ISF) compute misfdep 1201 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1202 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1203 END WHERE 1204 1205 ! Compute misfdep for ocean points (i.e. first wet level) 1206 ! find the first ocean level such that the first level thickness 1207 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1208 ! e3t_0 is the reference level thickness 1209 DO jk = 2, jpkm1 1210 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1211 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1212 END DO 1213 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1214 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1215 END WHERE 1216 1217 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1218 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1219 misfdep = 0; risfdep = 0.0_wp; 1220 mbathy = 0; bathy = 0.0_wp; 1221 END WHERE 1222 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1223 misfdep = 0; risfdep = 0.0_wp; 1224 mbathy = 0; bathy = 0.0_wp; 1225 END WHERE 1226 1227 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1228 icompt = 0 1229 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1230 DO jl = 1, 10 1231 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1232 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1233 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1234 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1235 END WHERE 1236 WHERE (mbathy(:,:) <= 0) 1237 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1238 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1239 END WHERE 1240 IF( lk_mpp ) THEN 1241 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1242 CALL lbc_lnk( zbathy, 'T', 1. ) 1243 misfdep(:,:) = INT( zbathy(:,:) ) 1244 1245 CALL lbc_lnk( risfdep,'T', 1. ) 1246 CALL lbc_lnk( bathy, 'T', 1. ) 1247 1248 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1249 CALL lbc_lnk( zbathy, 'T', 1. ) 1250 mbathy(:,:) = INT( zbathy(:,:) ) 1251 ENDIF 1252 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1253 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1254 misfdep(jpi,:) = misfdep( 2 ,:) 1255 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1256 mbathy(jpi,:) = mbathy( 2 ,:) 1257 ENDIF 1258 1259 ! split last cell if possible (only where water column is 2 cell or less) 1260 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1261 IF ( .NOT. ln_iscpl) THEN 1262 DO jk = jpkm1, 1, -1 1263 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1264 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1265 mbathy(:,:) = jk 1266 bathy(:,:) = zmax 1267 END WHERE 1268 END DO 1269 END IF 1270 1271 ! split top cell if possible (only where water column is 2 cell or less) 1272 DO jk = 2, jpkm1 1273 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1274 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1275 misfdep(:,:) = jk 1276 risfdep(:,:) = zmax 1277 END WHERE 1278 END DO 1279 1280 1281 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1282 DO jj = 1, jpj 1283 DO ji = 1, jpi 1284 ! find the minimum change option: 1285 ! test bathy 1286 IF (risfdep(ji,jj) > 1) THEN 1287 IF ( .NOT. ln_iscpl ) THEN 1288 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1289 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1290 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1291 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1292 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1293 IF (zbathydiff <= zrisfdepdiff) THEN 1294 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1295 mbathy(ji,jj)= mbathy(ji,jj) + 1 1296 ELSE 1297 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1298 misfdep(ji,jj) = misfdep(ji,jj) - 1 1299 END IF 1300 ENDIF 1301 ELSE 1302 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1303 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1304 misfdep(ji,jj) = misfdep(ji,jj) - 1 1305 END IF 1306 END IF 1307 END IF 1308 END DO 1309 END DO 1310 1311 ! At least 2 levels for water thickness at T, U, and V point. 1312 DO jj = 1, jpj 1313 DO ji = 1, jpi 1314 ! find the minimum change option: 1315 ! test bathy 1316 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1317 IF ( .NOT. ln_iscpl ) THEN 1318 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1319 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1320 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1321 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1322 IF (zbathydiff <= zrisfdepdiff) THEN 1323 mbathy(ji,jj) = mbathy(ji,jj) + 1 1324 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1325 ELSE 1326 misfdep(ji,jj)= misfdep(ji,jj) - 1 1327 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1328 END IF 1329 ELSE 1330 misfdep(ji,jj)= misfdep(ji,jj) - 1 1331 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1332 END IF 1333 ENDIF 1334 END DO 1335 END DO 1336 1337 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1338 DO jj = 1, jpjm1 1339 DO ji = 1, jpim1 1340 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1341 IF ( .NOT. ln_iscpl ) THEN 1342 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1343 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1344 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1345 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1346 IF (zbathydiff <= zrisfdepdiff) THEN 1347 mbathy(ji,jj) = mbathy(ji,jj) + 1 1348 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1349 ELSE 1350 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1351 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1352 END IF 1353 ELSE 1354 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1355 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1356 END IF 1357 ENDIF 1358 END DO 1359 END DO 1360 1361 IF( lk_mpp ) THEN 1362 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1363 CALL lbc_lnk( zbathy, 'T', 1. ) 1364 misfdep(:,:) = INT( zbathy(:,:) ) 1365 1366 CALL lbc_lnk( risfdep,'T', 1. ) 1367 CALL lbc_lnk( bathy, 'T', 1. ) 1368 1369 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1370 CALL lbc_lnk( zbathy, 'T', 1. ) 1371 mbathy(:,:) = INT( zbathy(:,:) ) 1372 ENDIF 1373 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1374 DO jj = 1, jpjm1 1375 DO ji = 1, jpim1 1376 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1377 IF ( .NOT. ln_iscpl ) THEN 1378 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1379 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1380 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1381 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1382 IF (zbathydiff <= zrisfdepdiff) THEN 1383 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1384 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1385 ELSE 1386 misfdep(ji,jj) = misfdep(ji,jj) - 1 1387 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1388 END IF 1389 ELSE 1390 misfdep(ji,jj) = misfdep(ji,jj) - 1 1391 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1392 END IF 1393 ENDIF 1394 END DO 1395 END DO 1396 1397 1398 IF( lk_mpp ) THEN 1399 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1400 CALL lbc_lnk( zbathy, 'T', 1. ) 1401 misfdep(:,:) = INT( zbathy(:,:) ) 1402 1403 CALL lbc_lnk( risfdep,'T', 1. ) 1404 CALL lbc_lnk( bathy, 'T', 1. ) 1405 1406 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1407 CALL lbc_lnk( zbathy, 'T', 1. ) 1408 mbathy(:,:) = INT( zbathy(:,:) ) 1409 ENDIF 1410 1411 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1412 DO jj = 1, jpjm1 1413 DO ji = 1, jpim1 1414 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1415 IF ( .NOT. ln_iscpl ) THEN 1416 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1417 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1418 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1419 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1420 IF (zbathydiff <= zrisfdepdiff) THEN 1421 mbathy(ji,jj) = mbathy(ji,jj) + 1 1422 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1423 ELSE 1424 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1425 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1426 END IF 1427 ELSE 1428 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1429 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1430 ENDIF 1431 ENDIF 1432 ENDDO 1433 ENDDO 1434 1435 IF( lk_mpp ) THEN 1436 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1437 CALL lbc_lnk( zbathy, 'T', 1. ) 1438 misfdep(:,:) = INT( zbathy(:,:) ) 1439 1440 CALL lbc_lnk( risfdep,'T', 1. ) 1441 CALL lbc_lnk( bathy, 'T', 1. ) 1442 1443 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1444 CALL lbc_lnk( zbathy, 'T', 1. ) 1445 mbathy(:,:) = INT( zbathy(:,:) ) 1446 ENDIF 1447 1448 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1449 DO jj = 1, jpjm1 1450 DO ji = 1, jpim1 1451 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1452 IF ( .NOT. ln_iscpl ) THEN 1453 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1454 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1455 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1456 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1457 IF (zbathydiff <= zrisfdepdiff) THEN 1458 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1459 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1460 ELSE 1461 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1462 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1463 END IF 1464 ELSE 1465 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1466 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1467 ENDIF 1468 ENDIF 1469 ENDDO 1470 ENDDO 1471 1472 IF( lk_mpp ) THEN 1473 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1474 CALL lbc_lnk( zbathy, 'T', 1. ) 1475 misfdep(:,:) = INT( zbathy(:,:) ) 1476 1477 CALL lbc_lnk( risfdep,'T', 1. ) 1478 CALL lbc_lnk( bathy, 'T', 1. ) 1479 1480 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1481 CALL lbc_lnk( zbathy, 'T', 1. ) 1482 mbathy(:,:) = INT( zbathy(:,:) ) 1483 ENDIF 1484 END DO 1485 ! end dig bathy/ice shelf to be compatible 1486 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1487 DO jl = 1,20 1488 1489 ! remove single point "bay" on isf coast line in the ice shelf draft' 1490 DO jk = 2, jpk 1491 WHERE (misfdep==0) misfdep=jpk 1492 zmask=0._wp 1493 WHERE (misfdep <= jk) zmask=1 1494 DO jj = 2, jpjm1 1495 DO ji = 2, jpim1 1496 IF (misfdep(ji,jj) == jk) THEN 1497 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1498 IF (ibtest <= 1) THEN 1499 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1500 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1501 END IF 1502 END IF 1503 END DO 1504 END DO 1505 END DO 1506 WHERE (misfdep==jpk) 1507 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1508 END WHERE 1509 IF( lk_mpp ) THEN 1510 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1511 CALL lbc_lnk( zbathy, 'T', 1. ) 1512 misfdep(:,:) = INT( zbathy(:,:) ) 1513 1514 CALL lbc_lnk( risfdep,'T', 1. ) 1515 CALL lbc_lnk( bathy, 'T', 1. ) 1516 1517 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1518 CALL lbc_lnk( zbathy, 'T', 1. ) 1519 mbathy(:,:) = INT( zbathy(:,:) ) 1520 ENDIF 1521 1522 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1523 DO jk = jpk,1,-1 1524 zmask=0._wp 1525 WHERE (mbathy >= jk ) zmask=1 1526 DO jj = 2, jpjm1 1527 DO ji = 2, jpim1 1528 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1529 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1530 IF (ibtest <= 1) THEN 1531 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1532 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1533 END IF 1534 END IF 1535 END DO 1536 END DO 1537 END DO 1538 WHERE (mbathy==0) 1539 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1540 END WHERE 1541 IF( lk_mpp ) THEN 1542 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1543 CALL lbc_lnk( zbathy, 'T', 1. ) 1544 misfdep(:,:) = INT( zbathy(:,:) ) 1545 1546 CALL lbc_lnk( risfdep,'T', 1. ) 1547 CALL lbc_lnk( bathy, 'T', 1. ) 1548 1549 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1550 CALL lbc_lnk( zbathy, 'T', 1. ) 1551 mbathy(:,:) = INT( zbathy(:,:) ) 1552 ENDIF 1553 1554 ! fill hole in ice shelf 1555 zmisfdep = misfdep 1556 zrisfdep = risfdep 1557 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1558 DO jj = 2, jpjm1 1559 DO ji = 2, jpim1 1560 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1561 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1562 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1563 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1564 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1565 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1566 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1567 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1568 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1569 END IF 1570 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1571 misfdep(ji,jj) = ibtest 1572 risfdep(ji,jj) = gdepw_1d(ibtest) 1573 ENDIF 1574 ENDDO 1575 ENDDO 1576 1577 IF( lk_mpp ) THEN 1578 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1579 CALL lbc_lnk( zbathy, 'T', 1. ) 1580 misfdep(:,:) = INT( zbathy(:,:) ) 1581 1582 CALL lbc_lnk( risfdep, 'T', 1. ) 1583 CALL lbc_lnk( bathy, 'T', 1. ) 1584 1585 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1586 CALL lbc_lnk( zbathy, 'T', 1. ) 1587 mbathy(:,:) = INT( zbathy(:,:) ) 1588 ENDIF 1589 ! 1590 !! fill hole in bathymetry 1591 zmbathy (:,:)=mbathy (:,:) 1592 DO jj = 2, jpjm1 1593 DO ji = 2, jpim1 1594 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1595 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1596 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1597 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1598 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1599 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1600 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1601 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1602 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1603 END IF 1604 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1605 mbathy(ji,jj) = ibtest 1606 bathy(ji,jj) = gdepw_1d(ibtest+1) 1607 ENDIF 1608 END DO 1609 END DO 1610 IF( lk_mpp ) THEN 1611 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1612 CALL lbc_lnk( zbathy, 'T', 1. ) 1613 misfdep(:,:) = INT( zbathy(:,:) ) 1614 1615 CALL lbc_lnk( risfdep, 'T', 1. ) 1616 CALL lbc_lnk( bathy, 'T', 1. ) 1617 1618 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1619 CALL lbc_lnk( zbathy, 'T', 1. ) 1620 mbathy(:,:) = INT( zbathy(:,:) ) 1621 ENDIF 1622 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1623 DO jj = 1, jpjm1 1624 DO ji = 1, jpim1 1625 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1626 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1627 END IF 1628 END DO 1629 END DO 1630 IF( lk_mpp ) THEN 1631 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1632 CALL lbc_lnk( zbathy, 'T', 1. ) 1633 misfdep(:,:) = INT( zbathy(:,:) ) 1634 1635 CALL lbc_lnk( risfdep, 'T', 1. ) 1636 CALL lbc_lnk( bathy, 'T', 1. ) 1637 1638 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1639 CALL lbc_lnk( zbathy, 'T', 1. ) 1640 mbathy(:,:) = INT( zbathy(:,:) ) 1641 ENDIF 1642 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1643 DO jj = 1, jpjm1 1644 DO ji = 1, jpim1 1645 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1646 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1647 END IF 1648 END DO 1649 END DO 1650 IF( lk_mpp ) THEN 1651 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1652 CALL lbc_lnk( zbathy, 'T', 1. ) 1653 misfdep(:,:) = INT( zbathy(:,:) ) 1654 1655 CALL lbc_lnk( risfdep,'T', 1. ) 1656 CALL lbc_lnk( bathy, 'T', 1. ) 1657 1658 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1659 CALL lbc_lnk( zbathy, 'T', 1. ) 1660 mbathy(:,:) = INT( zbathy(:,:) ) 1661 ENDIF 1662 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1663 DO jj = 1, jpjm1 1664 DO ji = 1, jpi 1665 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1666 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1667 END IF 1668 END DO 1669 END DO 1670 IF( lk_mpp ) THEN 1671 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1672 CALL lbc_lnk( zbathy, 'T', 1. ) 1673 misfdep(:,:) = INT( zbathy(:,:) ) 1674 1675 CALL lbc_lnk( risfdep,'T', 1. ) 1676 CALL lbc_lnk( bathy, 'T', 1. ) 1677 1678 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1679 CALL lbc_lnk( zbathy, 'T', 1. ) 1680 mbathy(:,:) = INT( zbathy(:,:) ) 1681 ENDIF 1682 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1683 DO jj = 1, jpjm1 1684 DO ji = 1, jpi 1685 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1686 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1687 END IF 1688 END DO 1689 END DO 1690 IF( lk_mpp ) THEN 1691 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1692 CALL lbc_lnk( zbathy, 'T', 1. ) 1693 misfdep(:,:) = INT( zbathy(:,:) ) 1694 1695 CALL lbc_lnk( risfdep,'T', 1. ) 1696 CALL lbc_lnk( bathy, 'T', 1. ) 1697 1698 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1699 CALL lbc_lnk( zbathy, 'T', 1. ) 1700 mbathy(:,:) = INT( zbathy(:,:) ) 1701 ENDIF 1702 ! if not compatible after all check, mask T 1703 DO jj = 1, jpj 1704 DO ji = 1, jpi 1705 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1706 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1707 END IF 1708 END DO 1709 END DO 1710 1711 WHERE (mbathy(:,:) == 1) 1712 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1713 END WHERE 1714 END DO 1715 ! end check compatibility ice shelf/bathy 1716 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1717 WHERE (risfdep(:,:) <= 10._wp) 1718 misfdep = 1; risfdep = 0.0_wp; 1719 END WHERE 1720 1721 IF( icompt == 0 ) THEN 1722 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1723 ELSE 1724 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1725 ENDIF 1726 1727 ! compute scale factor and depth at T- and W- points 1005 1728 DO jj = 1, jpj 1006 1729 DO ji = 1, jpi … … 1024 1747 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1025 1748 ENDIF 1749 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1026 1750 !gm Bug? check the gdepw_1d 1027 1751 ! ... on ik … … 1029 1753 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1030 1754 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1031 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1032 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1033 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1034 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1755 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1756 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1035 1757 ! ... on ik+1 1036 1758 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1037 1759 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1038 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1039 1760 ENDIF 1040 1761 ENDIF … … 1062 1783 END DO 1063 1784 ! 1064 IF ( ln_isfcav ) THEN1065 1785 ! (ISF) Definition of e3t, u, v, w for ISF case 1066 1067 1068 1069 1070 1071 1786 DO jj = 1, jpj 1787 DO ji = 1, jpi 1788 ik = misfdep(ji,jj) 1789 IF( ik > 1 ) THEN ! ice shelf point only 1790 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1791 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1072 1792 !gm Bug? check the gdepw_0 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 e3w_0 (ji,jj,ik ) =2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1085 1793 ! ... on ik 1794 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1795 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1796 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1797 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1798 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1799 1800 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1801 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1802 ENDIF 1803 ! ... on ik / ik-1 1804 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1805 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1086 1806 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1087 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1807 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1808 ENDIF 1809 END DO 1810 END DO 1811 1812 it = 0 1813 DO jj = 1, jpj 1814 DO ji = 1, jpi 1815 ik = misfdep(ji,jj) 1816 IF( ik > 1 ) THEN ! ice shelf point only 1817 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1818 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1819 ! test 1820 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1821 IF( zdiff <= 0. .AND. lwp ) THEN 1822 it = it + 1 1823 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1824 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1825 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1826 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1088 1827 ENDIF 1089 END DO1828 ENDIF 1090 1829 END DO 1091 !1092 it = 01093 DO jj = 1, jpj1094 DO ji = 1, jpi1095 ik = misfdep(ji,jj)1096 IF( ik > 1 ) THEN ! ice shelf point only1097 e3tp (ji,jj) = e3t_0(ji,jj,ik )1098 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1099 ! test1100 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1101 IF( zdiff <= 0. .AND. lwp ) THEN1102 it = it + 11103 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1104 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1105 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1106 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1107 ENDIF1108 ENDIF1109 END DO1110 END DO1111 END IF1112 ! END (ISF)1113 1114 ! Scale factors and depth at U-, V-, UW and VW-points1115 DO jk = 1, jpk ! initialisation to z-scale factors1116 e3u_0 (:,:,jk) = e3t_1d(jk)1117 e3v_0 (:,:,jk) = e3t_1d(jk)1118 e3uw_0(:,:,jk) = e3w_1d(jk)1119 e3vw_0(:,:,jk) = e3w_1d(jk)1120 END DO1121 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1122 DO jj = 1, jpjm11123 DO ji = 1, fs_jpim1 ! vector opt.1124 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1125 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1126 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1127 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1128 END DO1129 END DO1130 END DO1131 IF ( ln_isfcav ) THEN1132 ! (ISF) define e3uw (adapted for 2 cells in the water column)1133 DO jj = 2, jpjm11134 DO ji = 2, fs_jpim1 ! vector opt.1135 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1136 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1137 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1138 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1139 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1140 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1141 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1142 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1143 END DO1144 END DO1145 END IF1146 1147 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1148 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1149 !1150 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1151 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1152 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1153 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1154 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1155 END DO1156 1157 ! Scale factor at F-point1158 DO jk = 1, jpk ! initialisation to z-scale factors1159 e3f_0(:,:,jk) = e3t_1d(jk)1160 END DO1161 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1162 DO jj = 1, jpjm11163 DO ji = 1, fs_jpim1 ! vector opt.1164 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1165 END DO1166 END DO1167 END DO1168 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1169 !1170 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1171 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1172 END DO1173 !!gm bug ? : must be a do loop with mj0,mj11174 !1175 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21176 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1177 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1178 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1179 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1180 1181 ! Control of the sign1182 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1183 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1184 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1185 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1186 1187 ! Compute gdep3w_0 (vertical sum of e3w)1188 IF ( ln_isfcav ) THEN ! if cavity1189 WHERE (misfdep == 0) misfdep = 11190 DO jj = 1,jpj1191 DO ji = 1,jpi1192 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1193 DO jk = 2, misfdep(ji,jj)1194 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1195 END DO1196 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1197 DO jk = misfdep(ji,jj) + 1, jpk1198 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1199 END DO1200 END DO1201 END DO1202 ELSE ! no cavity1203 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1204 DO jk = 2, jpk1205 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1206 END DO1207 END IF1208 ! ! ================= !1209 IF(lwp .AND. ll_print) THEN ! Control print !1210 ! ! ================= !1211 DO jj = 1,jpj1212 DO ji = 1, jpi1213 ik = MAX( mbathy(ji,jj), 1 )1214 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1215 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1216 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1217 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1218 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1219 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1220 END DO1221 END DO1222 WRITE(numout,*)1223 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1224 WRITE(numout,*)1225 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1226 WRITE(numout,*)1227 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1228 WRITE(numout,*)1229 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1230 WRITE(numout,*)1231 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1232 WRITE(numout,*)1233 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1234 ENDIF1235 !1236 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1237 !1238 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1239 !1240 END SUBROUTINE zgr_zps1241 1242 SUBROUTINE zgr_isf1243 !!----------------------------------------------------------------------1244 !! *** ROUTINE zgr_isf ***1245 !!1246 !! ** Purpose : check the bathymetry in levels1247 !!1248 !! ** Method : THe water column have to contained at least 2 cells1249 !! Bathymetry and isfdraft are modified (dig/close) to respect1250 !! this criterion.1251 !!1252 !!1253 !! ** Action : - test compatibility between isfdraft and bathy1254 !! - bathy and isfdraft are modified1255 !!----------------------------------------------------------------------1256 !!1257 INTEGER :: ji, jj, jk, jl ! dummy loop indices1258 INTEGER :: ik, it ! temporary integers1259 INTEGER :: id, jd, nprocd1260 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1261 LOGICAL :: ll_print ! Allow control print for debugging1262 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1263 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1264 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1265 REAL(wp) :: zdiff ! temporary scalar1266 REAL(wp) :: zrefdep ! temporary scalar1267 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1268 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1269 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1270 !!---------------------------------------------------------------------1271 !1272 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1273 !1274 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)1275 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )1276 1277 1278 ! (ISF) compute misfdep1279 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11280 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1281 END WHERE1282 1283 ! Compute misfdep for ocean points (i.e. first wet level)1284 ! find the first ocean level such that the first level thickness1285 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1286 ! e3t_0 is the reference level thickness1287 DO jk = 2, jpkm11288 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1289 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11290 1830 END DO 1291 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1292 risfdep(:,:) = 0. ; misfdep(:,:) = 11293 END WHERE1294 1295 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation1296 icompt = 01297 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1298 DO jl = 1, 101299 WHERE (bathy(:,:) .EQ. risfdep(:,:) )1300 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1301 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1302 END WHERE1303 WHERE (mbathy(:,:) <= 0)1304 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1305 mbathy (:,:) = 0; bathy (:,:) = 0._wp1306 ENDWHERE1307 IF( lk_mpp ) THEN1308 zbathy(:,:) = FLOAT( misfdep(:,:) )1309 CALL lbc_lnk( zbathy, 'T', 1. )1310 misfdep(:,:) = INT( zbathy(:,:) )1311 CALL lbc_lnk( risfdep, 'T', 1. )1312 CALL lbc_lnk( bathy, 'T', 1. )1313 zbathy(:,:) = FLOAT( mbathy(:,:) )1314 CALL lbc_lnk( zbathy, 'T', 1. )1315 mbathy(:,:) = INT( zbathy(:,:) )1316 ENDIF1317 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1318 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1319 misfdep(jpi,:) = misfdep( 2 ,:)1320 ENDIF1321 1322 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1323 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1324 mbathy(jpi,:) = mbathy( 2 ,:)1325 ENDIF1326 1327 ! split last cell if possible (only where water column is 2 cell or less)1328 DO jk = jpkm1, 1, -11329 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1330 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1331 mbathy(:,:) = jk1332 bathy(:,:) = zmax1333 END WHERE1334 END DO1335 1336 ! split top cell if possible (only where water column is 2 cell or less)1337 DO jk = 2, jpkm11338 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1339 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1340 misfdep(:,:) = jk1341 risfdep(:,:) = zmax1342 END WHERE1343 END DO1344 1345 1346 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1347 DO jj = 1, jpj1348 DO ji = 1, jpi1349 ! find the minimum change option:1350 ! test bathy1351 IF (risfdep(ji,jj) .GT. 1) THEN1352 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1353 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1354 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1355 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1356 1357 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN1358 IF (zbathydiff .LE. zrisfdepdiff) THEN1359 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1360 mbathy(ji,jj)= mbathy(ji,jj) + 11361 ELSE1362 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1363 misfdep(ji,jj) = misfdep(ji,jj) - 11364 END IF1365 END IF1366 END IF1367 END DO1368 END DO1369 1370 ! At least 2 levels for water thickness at T, U, and V point.1371 DO jj = 1, jpj1372 DO ji = 1, jpi1373 ! find the minimum change option:1374 ! test bathy1375 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1376 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1377 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1378 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1379 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1380 IF (zbathydiff .LE. zrisfdepdiff) THEN1381 mbathy(ji,jj) = mbathy(ji,jj) + 11382 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1383 ELSE1384 misfdep(ji,jj)= misfdep(ji,jj) - 11385 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1386 END IF1387 ENDIF1388 END DO1389 END DO1390 1391 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1392 DO jj = 1, jpjm11393 DO ji = 1, jpim11394 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1395 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1396 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1397 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1398 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1399 IF (zbathydiff .LE. zrisfdepdiff) THEN1400 mbathy(ji,jj) = mbathy(ji,jj) + 11401 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1402 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1403 ELSE1404 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11405 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1406 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1407 END IF1408 ENDIF1409 END DO1410 END DO1411 1412 IF( lk_mpp ) THEN1413 zbathy(:,:) = FLOAT( misfdep(:,:) )1414 CALL lbc_lnk( zbathy, 'T', 1. )1415 misfdep(:,:) = INT( zbathy(:,:) )1416 CALL lbc_lnk( risfdep, 'T', 1. )1417 CALL lbc_lnk( bathy, 'T', 1. )1418 zbathy(:,:) = FLOAT( mbathy(:,:) )1419 CALL lbc_lnk( zbathy, 'T', 1. )1420 mbathy(:,:) = INT( zbathy(:,:) )1421 ENDIF1422 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1423 DO jj = 1, jpjm11424 DO ji = 1, jpim11425 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1426 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1427 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1428 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1429 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1430 IF (zbathydiff .LE. zrisfdepdiff) THEN1431 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11432 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1433 ELSE1434 misfdep(ji,jj) = misfdep(ji,jj) - 11435 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1436 END IF1437 ENDIF1438 END DO1439 END DO1440 1441 1442 IF( lk_mpp ) THEN1443 zbathy(:,:) = FLOAT( misfdep(:,:) )1444 CALL lbc_lnk( zbathy, 'T', 1. )1445 misfdep(:,:) = INT( zbathy(:,:) )1446 CALL lbc_lnk( risfdep, 'T', 1. )1447 CALL lbc_lnk( bathy, 'T', 1. )1448 zbathy(:,:) = FLOAT( mbathy(:,:) )1449 CALL lbc_lnk( zbathy, 'T', 1. )1450 mbathy(:,:) = INT( zbathy(:,:) )1451 ENDIF1452 1453 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1454 DO jj = 1, jpjm11455 DO ji = 1, jpim11456 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1457 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1458 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1459 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1460 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1461 IF (zbathydiff .LE. zrisfdepdiff) THEN1462 mbathy(ji,jj) = mbathy(ji,jj) + 11463 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1464 ELSE1465 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11466 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1467 END IF1468 ENDIF1469 ENDDO1470 ENDDO1471 1472 IF( lk_mpp ) THEN1473 zbathy(:,:) = FLOAT( misfdep(:,:) )1474 CALL lbc_lnk( zbathy, 'T', 1. )1475 misfdep(:,:) = INT( zbathy(:,:) )1476 CALL lbc_lnk( risfdep, 'T', 1. )1477 CALL lbc_lnk( bathy, 'T', 1. )1478 zbathy(:,:) = FLOAT( mbathy(:,:) )1479 CALL lbc_lnk( zbathy, 'T', 1. )1480 mbathy(:,:) = INT( zbathy(:,:) )1481 ENDIF1482 1483 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1484 DO jj = 1, jpjm11485 DO ji = 1, jpim11486 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1487 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1488 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1489 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1490 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1491 IF (zbathydiff .LE. zrisfdepdiff) THEN1492 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11493 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1494 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1495 ELSE1496 misfdep(ji,jj) = misfdep(ji ,jj) - 11497 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1498 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1499 END IF1500 ENDIF1501 ENDDO1502 ENDDO1503 1504 IF( lk_mpp ) THEN1505 zbathy(:,:) = FLOAT( misfdep(:,:) )1506 CALL lbc_lnk( zbathy, 'T', 1. )1507 misfdep(:,:) = INT( zbathy(:,:) )1508 CALL lbc_lnk( risfdep, 'T', 1. )1509 CALL lbc_lnk( bathy, 'T', 1. )1510 zbathy(:,:) = FLOAT( mbathy(:,:) )1511 CALL lbc_lnk( zbathy, 'T', 1. )1512 mbathy(:,:) = INT( zbathy(:,:) )1513 ENDIF1514 END DO1515 ! end dig bathy/ice shelf to be compatible1516 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1517 DO jl = 1,201518 1519 ! remove single point "bay" on isf coast line in the ice shelf draft'1520 DO jk = 2, jpk1521 WHERE (misfdep==0) misfdep=jpk1522 zmask=01523 WHERE (misfdep .LE. jk) zmask=11524 DO jj = 2, jpjm11525 DO ji = 2, jpim11526 IF (misfdep(ji,jj) .EQ. jk) THEN1527 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1528 IF (ibtest .LE. 1) THEN1529 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11530 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1531 END IF1532 END IF1533 END DO1534 END DO1535 END DO1536 WHERE (misfdep==jpk)1537 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1538 END WHERE1539 IF( lk_mpp ) THEN1540 zbathy(:,:) = FLOAT( misfdep(:,:) )1541 CALL lbc_lnk( zbathy, 'T', 1. )1542 misfdep(:,:) = INT( zbathy(:,:) )1543 CALL lbc_lnk( risfdep, 'T', 1. )1544 CALL lbc_lnk( bathy, 'T', 1. )1545 zbathy(:,:) = FLOAT( mbathy(:,:) )1546 CALL lbc_lnk( zbathy, 'T', 1. )1547 mbathy(:,:) = INT( zbathy(:,:) )1548 ENDIF1549 1550 ! remove single point "bay" on bathy coast line beneath an ice shelf'1551 DO jk = jpk,1,-11552 zmask=01553 WHERE (mbathy .GE. jk ) zmask=11554 DO jj = 2, jpjm11555 DO ji = 2, jpim11556 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1557 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1558 IF (ibtest .LE. 1) THEN1559 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11560 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01561 END IF1562 END IF1563 END DO1564 END DO1565 END DO1566 WHERE (mbathy==0)1567 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1568 END WHERE1569 IF( lk_mpp ) THEN1570 zbathy(:,:) = FLOAT( misfdep(:,:) )1571 CALL lbc_lnk( zbathy, 'T', 1. )1572 misfdep(:,:) = INT( zbathy(:,:) )1573 CALL lbc_lnk( risfdep, 'T', 1. )1574 CALL lbc_lnk( bathy, 'T', 1. )1575 zbathy(:,:) = FLOAT( mbathy(:,:) )1576 CALL lbc_lnk( zbathy, 'T', 1. )1577 mbathy(:,:) = INT( zbathy(:,:) )1578 ENDIF1579 1580 ! fill hole in ice shelf1581 zmisfdep = misfdep1582 zrisfdep = risfdep1583 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1584 DO jj = 2, jpjm11585 DO ji = 2, jpim11586 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1587 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1588 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1589 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1590 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1591 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1592 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1593 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1594 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1595 END IF1596 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1597 misfdep(ji,jj) = ibtest1598 risfdep(ji,jj) = gdepw_1d(ibtest)1599 ENDIF1600 ENDDO1601 ENDDO1602 1603 IF( lk_mpp ) THEN1604 zbathy(:,:) = FLOAT( misfdep(:,:) )1605 CALL lbc_lnk( zbathy, 'T', 1. )1606 misfdep(:,:) = INT( zbathy(:,:) )1607 CALL lbc_lnk( risfdep, 'T', 1. )1608 CALL lbc_lnk( bathy, 'T', 1. )1609 zbathy(:,:) = FLOAT( mbathy(:,:) )1610 CALL lbc_lnk( zbathy, 'T', 1. )1611 mbathy(:,:) = INT( zbathy(:,:) )1612 ENDIF1613 !1614 !! fill hole in bathymetry1615 zmbathy (:,:)=mbathy (:,:)1616 DO jj = 2, jpjm11617 DO ji = 2, jpim11618 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1619 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1620 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1621 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01622 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01623 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01624 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1625 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1626 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1627 END IF1628 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1629 mbathy(ji,jj) = ibtest1630 bathy(ji,jj) = gdepw_1d(ibtest+1)1631 ENDIF1632 END DO1633 END DO1634 IF( lk_mpp ) THEN1635 zbathy(:,:) = FLOAT( misfdep(:,:) )1636 CALL lbc_lnk( zbathy, 'T', 1. )1637 misfdep(:,:) = INT( zbathy(:,:) )1638 CALL lbc_lnk( risfdep, 'T', 1. )1639 CALL lbc_lnk( bathy, 'T', 1. )1640 zbathy(:,:) = FLOAT( mbathy(:,:) )1641 CALL lbc_lnk( zbathy, 'T', 1. )1642 mbathy(:,:) = INT( zbathy(:,:) )1643 ENDIF1644 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1645 DO jj = 1, jpjm11646 DO ji = 1, jpim11647 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1648 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1649 END IF1650 END DO1651 END DO1652 IF( lk_mpp ) THEN1653 zbathy(:,:) = FLOAT( misfdep(:,:) )1654 CALL lbc_lnk( zbathy, 'T', 1. )1655 misfdep(:,:) = INT( zbathy(:,:) )1656 CALL lbc_lnk( risfdep, 'T', 1. )1657 CALL lbc_lnk( bathy, 'T', 1. )1658 zbathy(:,:) = FLOAT( mbathy(:,:) )1659 CALL lbc_lnk( zbathy, 'T', 1. )1660 mbathy(:,:) = INT( zbathy(:,:) )1661 ENDIF1662 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1663 DO jj = 1, jpjm11664 DO ji = 1, jpim11665 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1666 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1667 END IF1668 END DO1669 END DO1670 IF( lk_mpp ) THEN1671 zbathy(:,:) = FLOAT( misfdep(:,:) )1672 CALL lbc_lnk( zbathy, 'T', 1. )1673 misfdep(:,:) = INT( zbathy(:,:) )1674 CALL lbc_lnk( risfdep, 'T', 1. )1675 CALL lbc_lnk( bathy, 'T', 1. )1676 zbathy(:,:) = FLOAT( mbathy(:,:) )1677 CALL lbc_lnk( zbathy, 'T', 1. )1678 mbathy(:,:) = INT( zbathy(:,:) )1679 ENDIF1680 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1681 DO jj = 1, jpjm11682 DO ji = 1, jpi1683 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1684 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1685 END IF1686 END DO1687 END DO1688 IF( lk_mpp ) THEN1689 zbathy(:,:) = FLOAT( misfdep(:,:) )1690 CALL lbc_lnk( zbathy, 'T', 1. )1691 misfdep(:,:) = INT( zbathy(:,:) )1692 CALL lbc_lnk( risfdep, 'T', 1. )1693 CALL lbc_lnk( bathy, 'T', 1. )1694 zbathy(:,:) = FLOAT( mbathy(:,:) )1695 CALL lbc_lnk( zbathy, 'T', 1. )1696 mbathy(:,:) = INT( zbathy(:,:) )1697 ENDIF1698 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1699 DO jj = 1, jpjm11700 DO ji = 1, jpi1701 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1702 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1703 END IF1704 END DO1705 END DO1706 IF( lk_mpp ) THEN1707 zbathy(:,:) = FLOAT( misfdep(:,:) )1708 CALL lbc_lnk( zbathy, 'T', 1. )1709 misfdep(:,:) = INT( zbathy(:,:) )1710 CALL lbc_lnk( risfdep, 'T', 1. )1711 CALL lbc_lnk( bathy, 'T', 1. )1712 zbathy(:,:) = FLOAT( mbathy(:,:) )1713 CALL lbc_lnk( zbathy, 'T', 1. )1714 mbathy(:,:) = INT( zbathy(:,:) )1715 ENDIF1716 ! if not compatible after all check, mask T1717 DO jj = 1, jpj1718 DO ji = 1, jpi1719 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1720 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1721 END IF1722 END DO1723 END DO1724 1725 WHERE (mbathy(:,:) == 1)1726 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1727 END WHERE1728 END DO1729 ! end check compatibility ice shelf/bathy1730 ! remove very shallow ice shelf (less than ~ 10m if 75L)1731 WHERE (misfdep(:,:) <= 5)1732 misfdep = 1; risfdep = 0.0_wp;1733 END WHERE1734 1735 IF( icompt == 0 ) THEN1736 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1737 ELSE1738 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1739 ENDIF1740 1831 1741 1832 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1742 1833 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1743 1744 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1745 1746 END SUBROUTINE 1834 ! 1835 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1836 ! 1837 END SUBROUTINE zgr_isf 1838 1747 1839 1748 1840 SUBROUTINE zgr_sco … … 1790 1882 !! 1791 1883 !!---------------------------------------------------------------------- 1792 !1793 1884 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1794 1885 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1799 1890 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1800 1891 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1801 1892 !! 1802 1893 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1803 1894 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1804 1895 !!---------------------------------------------------------------------- 1805 1896 ! 1806 1897 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1807 1898 ! 1808 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1899 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1809 1900 ! 1810 1901 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1851 1942 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1852 1943 1853 DO jj = 1, jpj 1854 DO ji = 1, jpi 1855 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1856 END DO 1857 END DO 1944 IF( .NOT.ln_wd ) THEN 1945 DO jj = 1, jpj 1946 DO ji = 1, jpi 1947 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1948 END DO 1949 END DO 1950 END IF 1858 1951 ! ! ============================= 1859 1952 ! ! Define the envelop bathymetry (hbatt) … … 1862 1955 zenv(:,:) = bathy(:,:) 1863 1956 ! 1957 IF( .NOT.ln_wd ) THEN 1864 1958 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1865 DO jj = 1, jpj 1866 DO ji = 1, jpi 1867 IF( bathy(ji,jj) == 0._wp ) THEN 1868 iip1 = MIN( ji+1, jpi ) 1869 ijp1 = MIN( jj+1, jpj ) 1870 iim1 = MAX( ji-1, 1 ) 1871 ijm1 = MAX( jj-1, 1 ) 1872 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1873 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1874 zenv(ji,jj) = rn_sbot_min 1875 ENDIF 1876 ENDIF 1877 END DO 1878 END DO 1959 DO jj = 1, jpj 1960 DO ji = 1, jpi 1961 IF( bathy(ji,jj) == 0._wp ) THEN 1962 iip1 = MIN( ji+1, jpi ) 1963 ijp1 = MIN( jj+1, jpj ) 1964 iim1 = MAX( ji-1, 1 ) 1965 ijm1 = MAX( jj-1, 1 ) 1966 !!gm BUG fix see ticket #1617 1967 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1968 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1969 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) & 1970 & zenv(ji,jj) = rn_sbot_min 1971 !!gm 1972 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1973 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1974 !!gm zenv(ji,jj) = rn_sbot_min 1975 !!gm ENDIF 1976 !!gm end 1977 ENDIF 1978 END DO 1979 END DO 1980 END IF 1981 1879 1982 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1880 1983 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1964 2067 ENDIF 1965 2068 ! 1966 IF(lwp) THEN ! Control print1967 WRITE(numout,*)1968 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1969 WRITE(numout,*)1970 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1971 IF( nprint == 1 ) THEN1972 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1973 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1974 ENDIF1975 ENDIF1976 1977 2069 ! ! ============================== 1978 2070 ! ! hbatu, hbatv, hbatf fields … … 1980 2072 IF(lwp) THEN 1981 2073 WRITE(numout,*) 1982 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2074 IF( .NOT.ln_wd ) THEN 2075 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2076 ELSE 2077 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2078 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2079 ENDIF 1983 2080 ENDIF 1984 2081 hbatu(:,:) = rn_sbot_min … … 1993 2090 END DO 1994 2091 END DO 2092 2093 IF( ln_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points 2094 DO jj = 1, jpj 2095 DO ji = 1, jpi 2096 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2097 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2098 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2099 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2100 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2101 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2102 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2103 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2104 END DO 2105 END DO 2106 END IF 1995 2107 ! 1996 2108 ! Apply lateral boundary condition … … 2000 2112 DO ji = 1, jpi 2001 2113 IF( hbatu(ji,jj) == 0._wp ) THEN 2114 !No worries about the following line when ln_wd == .true. 2002 2115 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2003 2116 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 2025 2138 2026 2139 !!bug: key_helsinki a verifer 2027 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2028 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2029 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2030 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2140 IF( .NOT.ln_wd ) THEN 2141 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2142 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2143 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2144 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2145 END IF 2031 2146 2032 2147 IF( nprint == 1 .AND. lwp ) THEN … … 2069 2184 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2070 2185 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2071 2072 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2073 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2074 ! 2075 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2076 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2077 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2078 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2079 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2080 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2081 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2186 ! 2187 IF( .NOT.ln_wd ) THEN 2188 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2189 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2190 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2191 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2192 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2193 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2194 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2195 END IF 2082 2196 2083 2197 #if defined key_agrif 2084 ! Ensure meaningful vertical scale factors in ghost lines/columns 2085 IF( .NOT. Agrif_Root() ) THEN 2086 ! 2087 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2088 e3u_0(1,:,:) = e3u_0(2,:,:) 2089 ENDIF 2090 ! 2091 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2092 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2093 ENDIF 2094 ! 2095 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2096 e3v_0(:,1,:) = e3v_0(:,2,:) 2097 ENDIF 2098 ! 2099 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2100 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2101 ENDIF 2102 ! 2103 ENDIF 2198 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2199 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2200 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2201 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2202 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2203 ENDIF 2104 2204 #endif 2105 2205 2106 fsdept(:,:,:) = gdept_0 (:,:,:) 2107 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2108 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2109 fse3t (:,:,:) = e3t_0 (:,:,:) 2110 fse3u (:,:,:) = e3u_0 (:,:,:) 2111 fse3v (:,:,:) = e3v_0 (:,:,:) 2112 fse3f (:,:,:) = e3f_0 (:,:,:) 2113 fse3w (:,:,:) = e3w_0 (:,:,:) 2114 fse3uw(:,:,:) = e3uw_0 (:,:,:) 2115 fse3vw(:,:,:) = e3vw_0 (:,:,:) 2206 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2207 !!gm and only that !!!!! 2208 !!gm THIS should be removed from here ! 2209 gdept_n(:,:,:) = gdept_0(:,:,:) 2210 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2211 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2212 e3t_n (:,:,:) = e3t_0 (:,:,:) 2213 e3u_n (:,:,:) = e3u_0 (:,:,:) 2214 e3v_n (:,:,:) = e3v_0 (:,:,:) 2215 e3f_n (:,:,:) = e3f_0 (:,:,:) 2216 e3w_n (:,:,:) = e3w_0 (:,:,:) 2217 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2218 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2219 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2220 !! gm end 2116 2221 !! 2117 2222 ! HYBRID : … … 2119 2224 DO ji = 1, jpi 2120 2225 DO jk = 1, jpkm1 2121 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2122 END DO 2123 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2226 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2227 END DO 2228 IF( ln_wd ) THEN 2229 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2230 mbathy(ji,jj) = 0 2231 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2232 mbathy(ji,jj) = 2 2233 ENDIF 2234 ELSE 2235 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2236 ENDIF 2124 2237 END DO 2125 2238 END DO … … 2130 2243 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2131 2244 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2132 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2133 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2134 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2135 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2245 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2246 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2247 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2248 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2136 2249 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2137 2250 2138 2251 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2139 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2140 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2141 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2142 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2252 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2253 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2254 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2255 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2143 2256 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2144 2257 ENDIF … … 2172 2285 END DO 2173 2286 ENDIF 2174 2175 !================================================================================2176 ! check the coordinate makes sense2177 !================================================================================2287 ! 2288 !================================================================================ 2289 ! check the coordinate makes sense 2290 !================================================================================ 2178 2291 DO ji = 1, jpi 2179 2292 DO jj = 1, jpj 2180 2293 ! 2181 2294 IF( hbatt(ji,jj) > 0._wp) THEN 2182 2295 DO jk = 1, mbathy(ji,jj) 2183 2296 ! check coordinate is monotonically increasing 2184 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2297 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2185 2298 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2186 2299 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2187 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2188 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2300 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2301 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2189 2302 CALL ctl_stop( ctmp1 ) 2190 2303 ENDIF 2191 2304 ! and check it has never gone negative 2192 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2305 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2193 2306 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2194 2307 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2195 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2196 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2308 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2309 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2197 2310 CALL ctl_stop( ctmp1 ) 2198 2311 ENDIF 2199 2312 ! and check it never exceeds the total depth 2200 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2313 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2201 2314 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2202 2315 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2203 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2316 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2204 2317 CALL ctl_stop( ctmp1 ) 2205 2318 ENDIF 2206 2319 END DO 2207 2320 ! 2208 2321 DO jk = 1, mbathy(ji,jj)-1 2209 2322 ! and check it never exceeds the total depth 2210 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2323 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2211 2324 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2212 2325 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2213 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2326 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2214 2327 CALL ctl_stop( ctmp1 ) 2215 2328 ENDIF 2216 2329 END DO 2217 2218 2330 ENDIF 2219 2220 2331 END DO 2221 2332 END DO … … 2227 2338 END SUBROUTINE zgr_sco 2228 2339 2229 !!====================================================================== 2340 2230 2341 SUBROUTINE s_sh94() 2231 2232 2342 !!---------------------------------------------------------------------- 2233 2343 !! *** ROUTINE s_sh94 *** … … 2240 2350 !! Reference : Song and Haidvogel 1994. 2241 2351 !!---------------------------------------------------------------------- 2242 !2243 2352 INTEGER :: ji, jj, jk ! dummy loop argument 2244 2353 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2354 REAL(wp) :: ztmpu, ztmpv, ztmpf 2355 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2245 2356 ! 2246 2357 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2247 2358 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2248 2249 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2250 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2359 !!---------------------------------------------------------------------- 2360 2361 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2362 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2251 2363 2252 2364 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2254 2366 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2255 2367 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2256 2368 ! 2257 2369 DO ji = 1, jpi 2258 2370 DO jj = 1, jpj 2259 2371 ! 2260 2372 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2261 2373 DO jk = 1, jpk … … 2286 2398 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2287 2399 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2288 gdept_0 2289 gdepw_0 2290 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2400 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2401 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2402 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2291 2403 END DO 2292 2404 ! … … 2296 2408 DO ji = 1, jpim1 2297 2409 DO jj = 1, jpjm1 2410 ! extended for Wetting/Drying case 2411 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2412 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2413 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2414 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2415 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2416 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2417 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2298 2418 DO jk = 1, jpk 2299 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2300 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2301 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2302 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2303 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2304 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2305 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2306 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2307 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2308 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2309 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2419 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2420 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2421 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2422 ELSE 2423 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2424 & / ztmpu 2425 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2426 & / ztmpu 2427 END IF 2428 2429 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2430 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2431 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2432 ELSE 2433 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2434 & / ztmpv 2435 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2436 & / ztmpv 2437 END IF 2438 2439 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2440 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj ,jk) + z_esigt3(ji+1,jj ,jk) & 2441 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2442 ELSE 2443 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2444 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2445 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2446 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2447 END IF 2448 2310 2449 ! 2311 2450 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 2320 2459 END DO 2321 2460 END DO 2322 2323 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2324 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2325 2461 ! 2462 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2463 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2464 ! 2326 2465 END SUBROUTINE s_sh94 2327 2466 2467 2328 2468 SUBROUTINE s_sf12 2329 2330 2469 !!---------------------------------------------------------------------- 2331 2470 !! *** ROUTINE s_sf12 *** … … 2343 2482 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2344 2483 !!---------------------------------------------------------------------- 2345 !2346 2484 INTEGER :: ji, jj, jk ! dummy loop argument 2347 2485 REAL(wp) :: zsmth ! smoothing around critical depth 2348 2486 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2487 REAL(wp) :: ztmpu, ztmpv, ztmpf 2488 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2349 2489 ! 2350 2490 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2351 2491 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2352 2492 !!---------------------------------------------------------------------- 2353 2493 ! 2354 2494 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2414 2554 2415 2555 DO jk = 1, jpk 2416 gdept_0 2417 gdepw_0 2418 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2556 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2557 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2558 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2419 2559 END DO 2420 2560 … … 2425 2565 DO jj=1,jpj-1 2426 2566 2427 DO jk = 1, jpk 2428 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 2429 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2430 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 2431 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2432 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 2433 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 2434 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2435 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 2436 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2437 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 2438 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2567 ! extend to suit for Wetting/Drying case 2568 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2569 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2570 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2571 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2572 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2573 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2574 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2575 DO jk = 1, jpk 2576 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2577 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2578 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2579 ELSE 2580 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2581 & / ztmpu 2582 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2583 & / ztmpu 2584 END IF 2585 2586 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2587 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2588 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2589 ELSE 2590 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2591 & / ztmpv 2592 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2593 & / ztmpv 2594 END IF 2595 2596 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2597 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2598 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2599 ELSE 2600 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2601 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2602 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2603 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2604 END IF 2605 2606 ! Code prior to wetting and drying option (for reference) 2607 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2608 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2609 ! 2610 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2611 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2612 ! 2613 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2614 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2615 ! 2616 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2617 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2618 ! 2619 !z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2620 ! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2621 ! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2622 ! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2623 ! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2439 2624 2440 2625 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) … … 2443 2628 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2444 2629 ! 2445 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2630 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2446 2631 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2447 2632 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2456 2641 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2457 2642 ! 2458 ! ! ============= 2459 2460 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2461 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2462 2643 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2644 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2645 ! 2463 2646 END SUBROUTINE s_sf12 2464 2647 2648 2465 2649 SUBROUTINE s_tanh() 2466 2467 2650 !!---------------------------------------------------------------------- 2468 2651 !! *** ROUTINE s_tanh*** … … 2474 2657 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2475 2658 !!---------------------------------------------------------------------- 2476 2477 INTEGER :: ji, jj, jk ! dummy loop argument 2659 INTEGER :: ji, jj, jk ! dummy loop argument 2478 2660 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2479 2480 2661 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2481 2662 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2482 2483 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2484 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2663 !!---------------------------------------------------------------------- 2664 2665 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2666 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2485 2667 2486 2668 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2512 2694 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2513 2695 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2514 gdept_0 2515 gdepw_0 2516 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2696 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2697 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2698 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2517 2699 END DO 2518 2700 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2531 2713 END DO 2532 2714 END DO 2533 2534 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2535 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2536 2715 ! 2716 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2717 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2718 ! 2537 2719 END SUBROUTINE s_tanh 2720 2538 2721 2539 2722 FUNCTION fssig( pk ) RESULT( pf ) … … 2607 2790 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2608 2791 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2609 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2610 REAL(wp), INTENT(in ) :: pzs ! surface box depth2611 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2612 REAL(wp) :: za1,za2,za3 ! local variables2613 REAL(wp) :: zn1,zn2 ! local variables2614 REAL(wp) :: za,zb,zx ! local variables2615 integer :: jk2616 !!----------------------------------------------------------------------2617 ! 2618 2619 zn1 = 1. /(jpk-1.)2620 zn2 = 1. - zn12621 2792 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2793 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2794 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2795 ! 2796 INTEGER :: jk ! dummy loop index 2797 REAL(wp) :: za1,za2,za3 ! local scalar 2798 REAL(wp) :: zn1,zn2 ! - - 2799 REAL(wp) :: za,zb,zx ! - - 2800 !!---------------------------------------------------------------------- 2801 ! 2802 zn1 = 1._wp / REAL( jpkm1, wp ) 2803 zn2 = 1._wp - zn1 2804 ! 2622 2805 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2623 2806 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2624 2807 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2625 2808 ! 2626 2809 za = pzb - za3*(pzs-za1)-za2 2627 2810 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2628 2811 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2629 2812 zx = 1.0_wp-za/2.0_wp-zb 2630 2813 ! 2631 2814 DO jk = 1, jpk 2632 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2633 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2634 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2815 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2816 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2817 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2635 2818 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2636 ENDDO 2637 2819 END DO 2638 2820 ! 2639 2821 END FUNCTION fgamma
Note: See TracChangeset
for help on using the changeset viewer.