- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4687 r6225 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 … … 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 19 21 !!---------------------------------------------------------------------- 20 22 … … 35 37 USE oce ! ocean variables 36 38 USE dom_oce ! ocean domain 39 USE wet_dry ! wetting and drying 37 40 USE closea ! closed seas 38 41 USE c1d ! 1D vertical configuration 42 ! 39 43 USE in_out_manager ! I/O manager 40 44 USE iom ! I/O library … … 72 76 73 77 !! * Substitutions 74 # include "domzgr_substitute.h90"75 78 # include "vectopt_loop_substitute.h90" 76 79 !!---------------------------------------------------------------------- … … 101 104 INTEGER :: ios 102 105 ! 103 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 104 107 !!---------------------------------------------------------------------- 105 108 ! … … 119 122 WRITE(numout,*) 'dom_zgr : vertical coordinate' 120 123 WRITE(numout,*) '~~~~~~~' 121 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 122 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 123 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 124 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 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 125 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 126 133 127 134 ioptio = 0 ! Check Vertical coordinate options … … 145 152 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 146 153 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 154 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 147 155 ! 148 156 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain … … 152 160 ! 153 161 IF( nprint == 1 .AND. lwp ) THEN 154 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )162 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 155 163 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 156 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )157 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &158 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &159 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &160 & ' 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(:,:,:) ) 161 169 162 170 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 163 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )164 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &165 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &166 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&167 & ' 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(:,:,:) ) 168 176 ENDIF 169 177 ! … … 216 224 & ppsur == pp_to_be_computed ) THEN 217 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 218 231 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 219 232 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 220 233 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 234 #endif 221 235 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 222 236 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) … … 233 247 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 234 248 WRITE(numout,*) ' Total depth :', zhmax 249 #if defined key_agrif 250 WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1) 251 #else 235 252 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 253 #endif 236 254 ELSE 237 255 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN … … 257 275 ! Reference z-coordinate (depth - scale factor at T- and W-points) 258 276 ! ====================== 259 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 260 281 za1 = zhmax / FLOAT(jpk-1) 282 #endif 261 283 DO jk = 1, jpk 262 284 zw = FLOAT( jk ) … … 294 316 gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero 295 317 ENDIF 318 319 IF ( ln_isfcav ) THEN 320 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 321 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 322 DO jk = 1, jpkm1 323 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 324 END DO 325 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 326 327 DO jk = 2, jpk 328 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 329 END DO 330 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 331 END IF 296 332 297 333 !!gm BUG in s-coordinate this does not work! … … 348 384 !! - bathy : meter bathymetry (in meters) 349 385 !!---------------------------------------------------------------------- 350 INTEGER :: ji, jj, j l, jk ! dummy loop indices386 INTEGER :: ji, jj, jk ! dummy loop indices 351 387 INTEGER :: inum ! temporary logical unit 388 INTEGER :: ierror ! error flag 352 389 INTEGER :: ii_bump, ij_bump, ih ! bump center position 353 390 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 354 391 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 355 392 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 356 INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data357 REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data393 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 394 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 358 395 !!---------------------------------------------------------------------- 359 396 ! 360 397 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 361 !362 CALL wrk_alloc( jpidta, jpjdta, idta )363 CALL wrk_alloc( jpidta, jpjdta, zdta )364 398 ! 365 399 IF(lwp) WRITE(numout,*) … … 370 404 ! ! ================== ! 371 405 ! ! global domain level and meter bathymetry (idta,zdta) 406 ! 407 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 408 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 409 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 410 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 372 411 ! 373 412 IF( ntopo == 0 ) THEN ! flat basin … … 450 489 END DO 451 490 END DO 491 risfdep(:,:)=0.e0 492 misfdep(:,:)=1 493 ! 494 DEALLOCATE( idta, zdta ) 452 495 ! 453 496 ! ! ================ ! … … 463 506 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 464 507 ! ! ===================== 465 IF( nn_cla == 0 ) THEN 466 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 467 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 468 DO ji = mi0(ii0), mi1(ii1) 469 DO jj = mj0(ij0), mj1(ij1) 470 mbathy(ji,jj) = 15 471 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 472 514 END DO 473 IF(lwp) WRITE(numout,*)474 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0475 !476 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open477 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)478 DO ji = mi0(ii0), mi1(ii1)479 DO jj = mj0(ij0), mj1(ij1)480 mbathy(ji,jj) = 12481 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 482 524 END DO 483 IF(lwp) WRITE(numout,*)484 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0485 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 486 528 ! 487 529 ENDIF … … 490 532 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 491 533 CALL iom_open ( 'bathy_meter.nc', inum ) 492 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 534 IF ( ln_isfcav ) THEN 535 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 536 ELSE 537 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 538 END IF 493 539 CALL iom_close( inum ) 494 540 ! 541 risfdep(:,:)=0._wp 542 misfdep(:,:)=1 543 IF ( ln_isfcav ) THEN 544 CALL iom_open ( 'isf_draft_meter.nc', inum ) 545 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 546 CALL iom_close( inum ) 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 555 END IF 556 ! 495 557 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 496 ! 497 IF( nn_cla == 0 ) THEN 498 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 499 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 500 DO ji = mi0(ii0), mi1(ii1) 501 DO jj = mj0(ij0), mj1(ij1) 502 bathy(ji,jj) = 284._wp 503 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 504 564 END DO 505 IF(lwp) WRITE(numout,*)506 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0507 !508 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open509 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)510 DO ji = mi0(ii0), mi1(ii1)511 DO jj = mj0(ij0), mj1(ij1)512 bathy(ji,jj) = 137._wp513 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 514 574 END DO 515 IF(lwp) WRITE(numout,*)516 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0517 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 518 578 ! 519 579 ENDIF … … 539 599 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 540 600 ENDIF 541 !542 CALL wrk_dealloc( jpidta, jpjdta, idta )543 CALL wrk_dealloc( jpidta, jpjdta, zdta )544 601 ! 545 602 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') … … 622 679 !! - update bathy : meter bathymetry (in meters) 623 680 !!---------------------------------------------------------------------- 624 !!625 681 INTEGER :: ji, jj, jl ! dummy loop indices 626 682 INTEGER :: icompt, ibtest, ikmax ! temporary integers 627 683 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 628 629 684 !!---------------------------------------------------------------------- 630 685 ! … … 723 778 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 724 779 ENDIF 725 726 IF( lwp .AND. nprint == 1 ) THEN ! control print727 WRITE(numout,*)728 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '729 WRITE(numout,*) ' ------------------'730 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )731 WRITE(numout,*)732 ENDIF733 780 ! 734 781 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 751 798 !! (min value = 1 over land) 752 799 !!---------------------------------------------------------------------- 753 !!754 800 INTEGER :: ji, jj ! dummy loop indices 755 801 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 784 830 785 831 832 SUBROUTINE zgr_top_level 833 !!---------------------------------------------------------------------- 834 !! *** ROUTINE zgr_top_level *** 835 !! 836 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 837 !! 838 !! ** Method : computes from misfdep with a minimum value of 1 839 !! 840 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 841 !! ocean level at t-, u- & v-points 842 !! (min value = 1) 843 !!---------------------------------------------------------------------- 844 INTEGER :: ji, jj ! dummy loop indices 845 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 846 !!---------------------------------------------------------------------- 847 ! 848 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 849 ! 850 CALL wrk_alloc( jpi, jpj, zmik ) 851 ! 852 IF(lwp) WRITE(numout,*) 853 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 854 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 855 ! 856 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 857 ! ! top k-index of W-level (=mikt) 858 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 859 DO ji = 1, jpim1 860 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 861 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 862 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 863 END DO 864 END DO 865 866 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 867 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 868 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 869 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 870 ! 871 CALL wrk_dealloc( jpi, jpj, zmik ) 872 ! 873 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 874 ! 875 END SUBROUTINE zgr_top_level 876 877 786 878 SUBROUTINE zgr_zco 787 879 !!---------------------------------------------------------------------- 788 880 !! *** ROUTINE zgr_zco *** 789 881 !! 790 !! ** Purpose : define the z-coordinate system882 !! ** Purpose : define the reference z-coordinate system 791 883 !! 792 884 !! ** Method : set 3D coord. arrays to reference 1D array … … 798 890 ! 799 891 DO jk = 1, jpk 800 gdept_0 801 gdepw_0 802 gde p3w_0(:,:,jk) = gdepw_1d(jk)803 e3t_0 804 e3u_0 805 e3v_0 806 e3f_0 807 e3w_0 808 e3uw_0 809 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) 810 902 END DO 811 903 ! … … 820 912 !! 821 913 !! ** Purpose : the depth and vertical scale factor in partial step 822 !! z-coordinate case914 !! reference z-coordinate case 823 915 !! 824 916 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 860 952 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 861 953 !!---------------------------------------------------------------------- 862 !!863 954 INTEGER :: ji, jj, jk ! dummy loop indices 864 INTEGER :: ik, it ! temporary integers 865 LOGICAL :: ll_print ! Allow control print for debugging 955 INTEGER :: ik, it, ikb, ikt ! temporary integers 866 956 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 867 957 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 868 REAL(wp) :: zmax ! Maximum depth869 958 REAL(wp) :: zdiff ! temporary scalar 870 REAL(wp) :: z refdep! temporary scalar959 REAL(wp) :: zmax ! temporary scalar 871 960 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 872 961 !!--------------------------------------------------------------------- … … 874 963 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 875 964 ! 876 CALL wrk_alloc( jpi, jpj, jpk,zprt )965 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 877 966 ! 878 967 IF(lwp) WRITE(numout,*) … … 880 969 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 881 970 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 882 883 ll_print = .FALSE. ! Local variable for debugging884 885 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth886 WRITE(numout,*)887 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'888 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )889 ENDIF890 891 971 892 972 ! bathymetry in level (from bathy_meter) … … 914 994 e3w_0 (:,:,jk) = e3w_1d (jk) 915 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 916 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 917 1728 DO jj = 1, jpj 918 1729 DO ji = 1, jpi … … 936 1747 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 937 1748 ENDIF 1749 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 938 1750 !gm Bug? check the gdepw_1d 939 1751 ! ... on ik 940 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 941 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 942 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 943 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 944 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 945 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 946 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1752 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1753 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1754 & / ( 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) 947 1757 ! ... on ik+1 948 1758 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 949 1759 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 950 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)951 1760 ENDIF 952 1761 ENDIF … … 973 1782 END DO 974 1783 END DO 975 976 ! Scale factors and depth at U-, V-, UW and VW-points 977 DO jk = 1, jpk ! initialisation to z-scale factors 978 e3u_0 (:,:,jk) = e3t_1d(jk) 979 e3v_0 (:,:,jk) = e3t_1d(jk) 980 e3uw_0(:,:,jk) = e3w_1d(jk) 981 e3vw_0(:,:,jk) = e3w_1d(jk) 982 END DO 983 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 984 DO jj = 1, jpjm1 985 DO ji = 1, fs_jpim1 ! vector opt. 986 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 987 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 988 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 989 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 990 END DO 991 END DO 992 END DO 993 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 994 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 995 ! 996 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 997 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 998 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 999 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1000 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1001 END DO 1002 1003 ! Scale factor at F-point 1004 DO jk = 1, jpk ! initialisation to z-scale factors 1005 e3f_0(:,:,jk) = e3t_1d(jk) 1006 END DO 1007 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1008 DO jj = 1, jpjm1 1009 DO ji = 1, fs_jpim1 ! vector opt. 1010 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1011 END DO 1012 END DO 1013 END DO 1014 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1015 ! 1016 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1017 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1018 END DO 1019 !!gm bug ? : must be a do loop with mj0,mj1 1020 ! 1021 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1022 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1023 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1024 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1025 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1026 1027 ! Control of the sign 1028 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1029 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1030 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1031 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1032 1033 ! Compute gdep3w_0 (vertical sum of e3w) 1034 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1035 DO jk = 2, jpk 1036 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1037 END DO 1038 1039 ! ! ================= ! 1040 IF(lwp .AND. ll_print) THEN ! Control print ! 1041 ! ! ================= ! 1042 DO jj = 1,jpj 1043 DO ji = 1, jpi 1044 ik = MAX( mbathy(ji,jj), 1 ) 1045 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1046 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1047 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1048 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1049 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1050 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1051 END DO 1052 END DO 1053 WRITE(numout,*) 1054 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1055 WRITE(numout,*) 1056 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1057 WRITE(numout,*) 1058 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1059 WRITE(numout,*) 1060 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1061 WRITE(numout,*) 1062 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1063 WRITE(numout,*) 1064 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1065 ENDIF 1066 ! 1067 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1068 ! 1069 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1070 ! 1071 END SUBROUTINE zgr_zps 1784 ! 1785 ! (ISF) Definition of e3t, u, v, w for ISF case 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) 1792 !gm Bug? check the gdepw_0 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) 1806 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 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) 1827 ENDIF 1828 ENDIF 1829 END DO 1830 END DO 1831 1832 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1833 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1834 ! 1835 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1836 ! 1837 END SUBROUTINE zgr_isf 1838 1072 1839 1073 1840 SUBROUTINE zgr_sco … … 1115 1882 !! 1116 1883 !!---------------------------------------------------------------------- 1117 !1118 1884 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1119 1885 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1124 1890 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1125 1891 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1126 1892 !! 1127 1893 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1128 1894 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1129 1895 !!---------------------------------------------------------------------- 1130 1896 ! 1131 1897 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1132 1898 ! 1133 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 ) 1134 1900 ! 1135 1901 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1176 1942 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1177 1943 1178 DO jj = 1, jpj 1179 DO ji = 1, jpi 1180 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1181 END DO 1182 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 1183 1951 ! ! ============================= 1184 1952 ! ! Define the envelop bathymetry (hbatt) … … 1187 1955 zenv(:,:) = bathy(:,:) 1188 1956 ! 1957 IF( .NOT.ln_wd ) THEN 1189 1958 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1190 DO jj = 1, jpj 1191 DO ji = 1, jpi 1192 IF( bathy(ji,jj) == 0._wp ) THEN 1193 iip1 = MIN( ji+1, jpi ) 1194 ijp1 = MIN( jj+1, jpj ) 1195 iim1 = MAX( ji-1, 1 ) 1196 ijm1 = MAX( jj-1, 1 ) 1197 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1198 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1199 zenv(ji,jj) = rn_sbot_min 1200 ENDIF 1201 ENDIF 1202 END DO 1203 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 1204 1982 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1205 1983 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1289 2067 ENDIF 1290 2068 ! 1291 IF(lwp) THEN ! Control print1292 WRITE(numout,*)1293 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1294 WRITE(numout,*)1295 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1296 IF( nprint == 1 ) THEN1297 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1298 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1299 ENDIF1300 ENDIF1301 1302 2069 ! ! ============================== 1303 2070 ! ! hbatu, hbatv, hbatf fields … … 1305 2072 IF(lwp) THEN 1306 2073 WRITE(numout,*) 1307 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 1308 2080 ENDIF 1309 2081 hbatu(:,:) = rn_sbot_min … … 1318 2090 END DO 1319 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 1320 2107 ! 1321 2108 ! Apply lateral boundary condition … … 1325 2112 DO ji = 1, jpi 1326 2113 IF( hbatu(ji,jj) == 0._wp ) THEN 2114 !No worries about the following line when ln_wd == .true. 1327 2115 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1328 2116 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 1350 2138 1351 2139 !!bug: key_helsinki a verifer 1352 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1353 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1354 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1355 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 1356 2146 1357 2147 IF( nprint == 1 .AND. lwp ) THEN … … 1394 2184 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 1395 2185 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1396 1397 fsdepw(:,:,:) = gdepw_0 (:,:,:) 1398 fsde3w(:,:,:) = gdep3w_0(:,:,:) 1399 ! 1400 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 1401 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 1402 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 1403 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 1404 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 1405 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 1406 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 1407 2196 1408 2197 #if defined key_agrif 1409 ! Ensure meaningful vertical scale factors in ghost lines/columns 1410 IF( .NOT. Agrif_Root() ) THEN 1411 ! 1412 IF((nbondi == -1).OR.(nbondi == 2)) THEN 1413 e3u_0(1,:,:) = e3u_0(2,:,:) 1414 ENDIF 1415 ! 1416 IF((nbondi == 1).OR.(nbondi == 2)) THEN 1417 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 1418 ENDIF 1419 ! 1420 IF((nbondj == -1).OR.(nbondj == 2)) THEN 1421 e3v_0(:,1,:) = e3v_0(:,2,:) 1422 ENDIF 1423 ! 1424 IF((nbondj == 1).OR.(nbondj == 2)) THEN 1425 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 1426 ENDIF 1427 ! 1428 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 1429 2204 #endif 1430 2205 1431 fsdept(:,:,:) = gdept_0 (:,:,:) 1432 fsdepw(:,:,:) = gdepw_0 (:,:,:) 1433 fsde3w(:,:,:) = gdep3w_0(:,:,:) 1434 fse3t (:,:,:) = e3t_0 (:,:,:) 1435 fse3u (:,:,:) = e3u_0 (:,:,:) 1436 fse3v (:,:,:) = e3v_0 (:,:,:) 1437 fse3f (:,:,:) = e3f_0 (:,:,:) 1438 fse3w (:,:,:) = e3w_0 (:,:,:) 1439 fse3uw(:,:,:) = e3uw_0 (:,:,:) 1440 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 1441 2221 !! 1442 2222 ! HYBRID : … … 1444 2224 DO ji = 1, jpi 1445 2225 DO jk = 1, jpkm1 1446 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1447 END DO 1448 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 1449 2237 END DO 1450 2238 END DO … … 1455 2243 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1456 2244 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 1457 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )1458 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 1459 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 1460 & ' 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 (:,:,:) ), & 1461 2249 & ' w ', MINVAL( e3w_0 (:,:,:) ) 1462 2250 1463 2251 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 1464 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )1465 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 1466 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 1467 & ' 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 (:,:,:) ), & 1468 2256 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 1469 2257 ENDIF … … 1497 2285 END DO 1498 2286 ENDIF 1499 1500 !================================================================================1501 ! check the coordinate makes sense1502 !================================================================================2287 ! 2288 !================================================================================ 2289 ! check the coordinate makes sense 2290 !================================================================================ 1503 2291 DO ji = 1, jpi 1504 2292 DO jj = 1, jpj 1505 2293 ! 1506 2294 IF( hbatt(ji,jj) > 0._wp) THEN 1507 2295 DO jk = 1, mbathy(ji,jj) 1508 2296 ! check coordinate is monotonically increasing 1509 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 1510 2298 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1511 2299 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1512 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)1513 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2300 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2301 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 1514 2302 CALL ctl_stop( ctmp1 ) 1515 2303 ENDIF 1516 2304 ! and check it has never gone negative 1517 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 1518 2306 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1519 2307 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1520 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)1521 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2308 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2309 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 1522 2310 CALL ctl_stop( ctmp1 ) 1523 2311 ENDIF 1524 2312 ! and check it never exceeds the total depth 1525 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2313 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 1526 2314 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1527 2315 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1528 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2316 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 1529 2317 CALL ctl_stop( ctmp1 ) 1530 2318 ENDIF 1531 2319 END DO 1532 2320 ! 1533 2321 DO jk = 1, mbathy(ji,jj)-1 1534 2322 ! and check it never exceeds the total depth 1535 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2323 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 1536 2324 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1537 2325 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1538 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2326 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 1539 2327 CALL ctl_stop( ctmp1 ) 1540 2328 ENDIF 1541 2329 END DO 1542 1543 2330 ENDIF 1544 1545 2331 END DO 1546 2332 END DO … … 1552 2338 END SUBROUTINE zgr_sco 1553 2339 1554 !!====================================================================== 2340 1555 2341 SUBROUTINE s_sh94() 1556 1557 2342 !!---------------------------------------------------------------------- 1558 2343 !! *** ROUTINE s_sh94 *** … … 1565 2350 !! Reference : Song and Haidvogel 1994. 1566 2351 !!---------------------------------------------------------------------- 1567 !1568 2352 INTEGER :: ji, jj, jk ! dummy loop argument 1569 2353 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2354 REAL(wp) :: ztmpu, ztmpv, ztmpf 2355 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 1570 2356 ! 1571 2357 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1572 2358 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1573 1574 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1575 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 ) 1576 2363 1577 2364 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 1579 2366 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 1580 2367 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 1581 2368 ! 1582 2369 DO ji = 1, jpi 1583 2370 DO jj = 1, jpj 1584 2371 ! 1585 2372 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1586 2373 DO jk = 1, jpk … … 1611 2398 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1612 2399 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1613 gdept_0 1614 gdepw_0 1615 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 ) 1616 2403 END DO 1617 2404 ! … … 1621 2408 DO ji = 1, jpim1 1622 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)) 1623 2418 DO jk = 1, jpk 1624 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1625 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1626 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1627 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1628 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1629 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1630 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1631 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1632 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1633 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1634 & / ( 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 1635 2449 ! 1636 2450 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 1645 2459 END DO 1646 2460 END DO 1647 1648 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )1649 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )1650 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 ! 1651 2465 END SUBROUTINE s_sh94 1652 2466 2467 1653 2468 SUBROUTINE s_sf12 1654 1655 2469 !!---------------------------------------------------------------------- 1656 2470 !! *** ROUTINE s_sf12 *** … … 1668 2482 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 1669 2483 !!---------------------------------------------------------------------- 1670 !1671 2484 INTEGER :: ji, jj, jk ! dummy loop argument 1672 2485 REAL(wp) :: zsmth ! smoothing around critical depth 1673 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 1674 2489 ! 1675 2490 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1676 2491 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1677 2492 !!---------------------------------------------------------------------- 1678 2493 ! 1679 2494 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 1739 2554 1740 2555 DO jk = 1, jpk 1741 gdept_0 1742 gdepw_0 1743 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) 1744 2559 END DO 1745 2560 … … 1750 2565 DO jj=1,jpj-1 1751 2566 1752 DO jk = 1, jpk 1753 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 1754 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1755 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 1756 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1757 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 1758 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 1759 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1760 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 1761 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1762 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 1763 ( 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) ) 1764 2624 1765 2625 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) … … 1768 2628 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 1769 2629 ! 1770 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) 1771 2631 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 1772 2632 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 1781 2641 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 1782 2642 ! 1783 ! ! ============= 1784 1785 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1786 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1787 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 ! 1788 2646 END SUBROUTINE s_sf12 1789 2647 2648 1790 2649 SUBROUTINE s_tanh() 1791 1792 2650 !!---------------------------------------------------------------------- 1793 2651 !! *** ROUTINE s_tanh*** … … 1799 2657 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1800 2658 !!---------------------------------------------------------------------- 1801 1802 INTEGER :: ji, jj, jk ! dummy loop argument 2659 INTEGER :: ji, jj, jk ! dummy loop argument 1803 2660 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1804 1805 2661 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 1806 2662 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 1807 1808 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 1809 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 ) 1810 2667 1811 2668 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 1837 2694 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1838 2695 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1839 gdept_0 1840 gdepw_0 1841 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 ) 1842 2699 END DO 1843 2700 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 1856 2713 END DO 1857 2714 END DO 1858 1859 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)1860 CALL wrk_dealloc( jpk, z_esigt, z_esigw)1861 2715 ! 2716 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2717 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2718 ! 1862 2719 END SUBROUTINE s_tanh 2720 1863 2721 1864 2722 FUNCTION fssig( pk ) RESULT( pf ) … … 1932 2790 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 1933 2791 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 1934 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth1935 REAL(wp), INTENT(in ) :: pzs ! surface box depth1936 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter1937 REAL(wp) :: za1,za2,za3 ! local variables1938 REAL(wp) :: zn1,zn2 ! local variables1939 REAL(wp) :: za,zb,zx ! local variables1940 integer :: jk1941 !!----------------------------------------------------------------------1942 ! 1943 1944 zn1 = 1. /(jpk-1.)1945 zn2 = 1. - zn11946 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 ! 1947 2805 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 1948 2806 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 1949 2807 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 1950 2808 ! 1951 2809 za = pzb - za3*(pzs-za1)-za2 1952 2810 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 1953 2811 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 1954 2812 zx = 1.0_wp-za/2.0_wp-zb 1955 2813 ! 1956 2814 DO jk = 1, jpk 1957 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &1958 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &1959 &(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) ) 1960 2818 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 1961 ENDDO 1962 2819 END DO 1963 2820 ! 1964 2821 END FUNCTION fgamma
Note: See TracChangeset
for help on using the changeset viewer.
