Changeset 12101 for utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domzgr.F90
- Timestamp:
- 2019-12-06T18:45:39+01:00 (4 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domzgr.F90
r12100 r12101 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 capability e19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capability 20 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 21 21 !!---------------------------------------------------------------------- … … 35 35 !! fgamma : Siddorn and Furner 2012 stretching function 36 36 !!--------------------------------------------------------------------- 37 USE oce ! ocean variables38 37 USE dom_oce ! ocean domain 39 ! USE closea ! closed seas40 38 ! 41 39 USE in_out_manager ! I/O manager … … 43 41 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 44 42 USE lib_mpp ! distributed memory computing library 45 USE wrk_nemo ! Memory allocation 46 USE timing ! Timing 43 USE lib_fortran 47 44 USE dombat 45 USE domisf 48 46 49 47 IMPLICIT NONE … … 51 49 52 50 PUBLIC dom_zgr ! called by dom_init.F90 53 51 ! 54 52 ! !!* Namelist namzgr_sco * 55 53 LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) … … 60 58 REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1) 61 59 REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates 62 INTEGER , PUBLIC :: ntopo !: = 0/1 ,compute/read the bathymetry file63 REAL(wp), PUBLIC :: e3zps_min !: miminum thickness for partial steps (meters)64 REAL(wp), PUBLIC :: e3zps_rat !: minimum thickness ration for partial steps65 INTEGER, PUBLIC :: nperio !: type of lateral boundary condition66 60 67 61 ! Song and Haidvogel 1994 stretching parameters … … 121 115 !!---------------------------------------------------------------------- 122 116 ! 123 ! IF( nn_timing == 1 ) CALL timing_start('dom_zgr')124 117 ! 125 118 REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate … … 152 145 IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 153 146 ! 147 IF ( ln_isfcav ) CALL zgr_isf_nam 154 148 ioptio = 0 155 149 IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1 … … 164 158 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 165 159 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 160 CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 166 161 ! 167 162 ! final adjustment of mbathy & check 168 163 ! ----------------------------------- 169 IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain170 CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points171 164 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 172 165 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points … … 175 168 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 176 169 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 177 & ' w ', MINVAL( gdepw_0(:,:,:) ) , '3w ', MINVAL( gde3w_0(:,:,:) )170 & ' w ', MINVAL( gdepw_0(:,:,:) ) 178 171 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 179 172 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & … … 182 175 183 176 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 184 & ' w ', MAXVAL( gdepw_0(:,:,:) ) , '3w ', MAXVAL( gde3w_0(:,:,:) )177 & ' w ', MAXVAL( gdepw_0(:,:,:) ) 185 178 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 186 179 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & … … 188 181 & ' w ', MAXVAL( e3w_0(:,:,:) ) 189 182 ENDIF 190 !191 ! IF( nn_timing == 1 ) CALL timing_stop('dom_zgr')192 183 ! 193 184 END SUBROUTINE dom_zgr … … 222 213 REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters 223 214 !!---------------------------------------------------------------------- 224 !225 ! IF( nn_timing == 1 ) CALL timing_start('zgr_z')226 215 ! 227 216 ! Set variables from parameters … … 322 311 IF ( ln_isfcav .OR. ln_e3_dep ) THEN ! e3. = dk[gdep] 323 312 ! 324 !==>>> need to be like this to compute the pressure gradient with ISF.325 ! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)326 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively327 !328 313 DO jk = 1, jpkm1 329 314 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) … … 354 339 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 355 340 END DO 356 !357 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_z')358 341 ! 359 342 END SUBROUTINE zgr_z … … 401 384 !!---------------------------------------------------------------------- 402 385 ! 403 ! IF( nn_timing == 1 ) CALL timing_start('zgr_bat')404 !405 386 IF(lwp) WRITE(numout,*) 406 387 IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry' … … 411 392 ! ! global domain level and meter bathymetry (idta,zdta) 412 393 ! 413 ALLOCATE( idta(jpi dta,jpjdta), STAT=ierror )394 ALLOCATE( idta(jpiglo,jpjglo), STAT=ierror ) 414 395 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 415 ALLOCATE( zdta(jpi dta,jpjdta), STAT=ierror )396 ALLOCATE( zdta(jpiglo,jpjglo), STAT=ierror ) 416 397 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 417 398 ! … … 439 420 IF(lwp) WRITE(numout,*) 440 421 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump' 441 ii_bump = jpi dta/ 2 ! i-index of the bump center442 ij_bump = jpj dta/ 2 ! j-index of the bump center422 ii_bump = jpiglo / 2 ! i-index of the bump center 423 ij_bump = jpjglo / 2 ! j-index of the bump center 443 424 r_bump = 50000._wp ! bump radius (meters) 444 425 h_bump = 2700._wp ! bump height (meters) … … 450 431 IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters' 451 432 ! 452 DO jj = 1, jpj dta! zdta :453 DO ji = 1, jpi dta433 DO jj = 1, jpjglo ! zdta : 434 DO ji = 1, jpiglo 454 435 zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 455 436 zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump … … 467 448 ENDIF 468 449 ENDIF 450 ! 469 451 ! ! set GLOBAL boundary conditions 470 ! ! Caution : idta on the global domain: use of jperio, not nperio471 452 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 472 453 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp 473 idta( : ,jpj dta) = 0 ; zdta( : ,jpjdta) = 0._wp454 idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp 474 455 ELSEIF( jperio == 2 ) THEN 475 456 idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 ) 476 idta( : ,jpj dta) = 0 ; zdta( : ,jpjdta) = 0._wp457 idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp 477 458 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp 478 idta(jpi dta, : ) = 0 ; zdta(jpidta, : ) = 0._wp459 idta(jpiglo, : ) = 0 ; zdta(jpiglo, : ) = 0._wp 479 460 ELSE 480 461 ih = 0 ; zh = 0._wp 481 462 IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce 482 463 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh 483 idta( : ,jpj dta) = ih ; zdta( : ,jpjdta) = zh464 idta( : ,jpjglo) = ih ; zdta( : ,jpjglo) = zh 484 465 idta( 1 , : ) = ih ; zdta( 1 , : ) = zh 485 idta(jpi dta, : ) = ih ; zdta(jpidta, : ) = zh466 idta(jpiglo, : ) = ih ; zdta(jpiglo, : ) = zh 486 467 ENDIF 487 468 … … 497 478 risfdep(:,:)=0.e0 498 479 misfdep(:,:)=1 499 !500 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code501 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN502 risfdep(:,:)=200.e0503 misfdep(:,:)=1504 ij0 = 1 ; ij1 = 40505 DO jj = mj0(ij0), mj1(ij1)506 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp507 END DO508 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp509 !510 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN511 !512 risfdep(:,:)=0.e0513 misfdep(:,:)=1514 ij0 = 1 ; ij1 = 40515 DO jj = mj0(ij0), mj1(ij1)516 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp517 END DO518 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp519 END IF520 480 ! 521 481 DEALLOCATE( idta, zdta ) … … 591 551 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 592 552 CALL iom_close( inum ) 593 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp594 595 ! set grounded point to 0596 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction)597 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin )598 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp599 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp600 END WHERE601 553 END IF 602 554 ! … … 633 585 ENDIF 634 586 ! 635 ! IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!636 !637 587 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 638 588 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 639 589 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth 640 590 ENDIF 641 zhmin = gdepw_1d(ik+1) 591 zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels 642 592 WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands 643 ELSE WHERE 593 ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans 644 594 END WHERE 645 595 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 646 596 ENDIF 647 597 ! 648 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')649 !650 598 END SUBROUTINE zgr_bat 651 652 653 SUBROUTINE zgr_bat_zoom654 !!----------------------------------------------------------------------655 !! *** ROUTINE zgr_bat_zoom ***656 !!657 !! ** Purpose : - Close zoom domain boundary if necessary658 !! - Suppress Med Sea from ORCA R2 and R05 arctic zoom659 !!660 !! ** Method :661 !!662 !! ** Action : - update mbathy: level bathymetry (in level index)663 !!----------------------------------------------------------------------664 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers665 !!----------------------------------------------------------------------666 !667 IF(lwp) WRITE(numout,*)668 IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'669 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'670 !671 ! Zoom domain672 ! ===========673 !674 ! Forced closed boundary if required675 IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0676 IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0677 IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0678 IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0679 !680 ! Configuration specific domain modifications681 ! (here, ORCA arctic configuration: suppress Med Sea)682 IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN683 SELECT CASE ( jp_cfg )684 ! ! =======================685 CASE ( 2 ) ! ORCA_R2 configuration686 ! ! =======================687 IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'688 ii0 = 141 ; ii1 = 162 ! Sea box i,j indices689 ij0 = 98 ; ij1 = 110690 ! ! =======================691 CASE ( 05 ) ! ORCA_R05 configuration692 ! ! =======================693 IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'694 ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe695 ij0 = 314 ; ij1 = 370696 END SELECT697 !698 mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe699 !700 ENDIF701 !702 END SUBROUTINE zgr_bat_zoom703 704 599 705 600 SUBROUTINE zgr_bat_ctl … … 727 622 INTEGER :: ji, jj, jl ! dummy loop indices 728 623 INTEGER :: icompt, ibtest, ikmax ! temporary integers 729 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 730 !!---------------------------------------------------------------------- 731 ! 732 ! IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl') 733 ! 734 CALL wrk_alloc( jpi, jpj, zbathy ) 624 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy 625 !!---------------------------------------------------------------------- 626 ! 627 ALLOCATE(zbathy(jpi,jpj)) 735 628 ! 736 629 IF(lwp) WRITE(numout,*) … … 743 636 icompt = 0 744 637 DO jl = 1, 2 745 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN638 IF( l_Iperio ) THEN 746 639 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 747 640 mbathy(jpi,:) = mbathy( 2 ,:) 748 641 ENDIF 642 zbathy(:,:) = FLOAT( mbathy(:,:) ) 643 CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 644 mbathy(:,:) = INT( zbathy(:,:) ) 645 749 646 DO jj = 2, jpjm1 750 647 DO ji = 2, jpim1 … … 760 657 END DO 761 658 END DO 762 ! IF( lk_mpp ) CALL mpp_sum( icompt ) 659 660 IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt ) 763 661 IF( icompt == 0 ) THEN 764 662 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' … … 766 664 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 767 665 ENDIF 768 IF( lk_mpp ) THEN 769 770 CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )771 772 ENDIF 666 667 zbathy(:,:) = FLOAT( mbathy(:,:) ) 668 CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 669 mbathy(:,:) = INT( zbathy(:,:) ) 670 773 671 ! ! East-west cyclic boundary conditions 774 IF( nperio == 0 ) THEN775 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio672 IF( jperio == 0 ) THEN 673 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio 776 674 IF( lk_mpp ) THEN 777 675 IF( nbondi == -1 .OR. nbondi == 2 ) THEN … … 790 688 ENDIF 791 689 ENDIF 792 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN793 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio690 ELSEIF( l_Iperio ) THEN 691 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio 794 692 mbathy( 1 ,:) = mbathy(jpim1,:) 795 693 mbathy(jpi,:) = mbathy( 2 ,:) 796 ELSEIF( nperio == 2 ) THEN797 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio694 ELSEIF( jperio == 2 ) THEN 695 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio 798 696 ELSE 799 697 IF(lwp) WRITE(numout,*) ' e r r o r' 800 IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio698 IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio 801 699 ! STOP 'dom_mba' 802 700 ENDIF 701 803 702 ! Boundary condition on mbathy 804 703 IF( .NOT.lk_mpp ) THEN … … 806 705 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 807 706 zbathy(:,:) = FLOAT( mbathy(:,:) ) 808 CALL lbc_lnk( ' toto',zbathy, 'T', 1._wp )707 CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 809 708 mbathy(:,:) = INT( zbathy(:,:) ) 810 709 ENDIF 710 811 711 ! Number of ocean level inferior or equal to jpkm1 812 ikmax = 0 813 DO jj = 1, jpj 814 DO ji = 1, jpi 815 ikmax = MAX( ikmax, mbathy(ji,jj) ) 816 END DO 817 END DO 818 !!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ??? 712 zbathy(:,:) = FLOAT( mbathy(:,:) ) 713 ikmax = glob_max( 'domzgr', zbathy(:,:) ) 714 819 715 IF( ikmax > jpkm1 ) THEN 820 716 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1' … … 825 721 ENDIF 826 722 ! 827 CALL wrk_dealloc( jpi, jpj, zbathy ) 828 ! 829 !! IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl') 723 DEALLOCATE( zbathy ) 830 724 ! 831 725 END SUBROUTINE zgr_bat_ctl … … 845 739 !!---------------------------------------------------------------------- 846 740 INTEGER :: ji, jj ! dummy loop indices 847 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 848 !!---------------------------------------------------------------------- 849 ! 850 ! IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level') 851 ! 852 CALL wrk_alloc( jpi, jpj, zmbk ) 741 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk 742 !!---------------------------------------------------------------------- 743 ! 744 ALLOCATE( zmbk(jpi,jpj) ) 853 745 ! 854 746 IF(lwp) WRITE(numout,*) … … 866 758 END DO 867 759 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 868 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('toto',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 869 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('toto',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 870 ! 871 CALL wrk_dealloc( jpi, jpj, zmbk ) 872 ! 873 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level') 760 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 761 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 762 ! 763 DEALLOCATE( zmbk ) 874 764 ! 875 765 END SUBROUTINE zgr_bot_level … … 889 779 !!---------------------------------------------------------------------- 890 780 INTEGER :: ji, jj ! dummy loop indices 891 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 892 !!---------------------------------------------------------------------- 893 ! 894 ! IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 895 ! 896 CALL wrk_alloc( jpi, jpj, zmik ) 781 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik 782 !!---------------------------------------------------------------------- 783 ! 784 ALLOCATE( zmik(jpi,jpj) ) 897 785 ! 898 786 IF(lwp) WRITE(numout,*) … … 911 799 912 800 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 913 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('toto',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 914 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('toto',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 915 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('toto',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 916 ! 917 CALL wrk_dealloc( jpi, jpj, zmik ) 918 ! 919 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 801 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 802 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 803 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 804 ! 805 DEALLOCATE( zmik ) 920 806 ! 921 807 END SUBROUTINE zgr_top_level … … 932 818 INTEGER :: jk 933 819 !!---------------------------------------------------------------------- 934 !935 ! IF( nn_timing == 1 ) CALL timing_start('zgr_zco')936 820 ! 937 821 DO jk = 1, jpk 938 822 gdept_0(:,:,jk) = gdept_1d(jk) 939 823 gdepw_0(:,:,jk) = gdepw_1d(jk) 940 gde3w_0(:,:,jk) = gdepw_1d(jk)941 824 e3t_0 (:,:,jk) = e3t_1d (jk) 942 825 e3u_0 (:,:,jk) = e3t_1d (jk) … … 947 830 e3vw_0 (:,:,jk) = e3w_1d (jk) 948 831 END DO 949 !950 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')951 832 ! 952 833 END SUBROUTINE zgr_zco … … 1004 885 REAL(wp) :: zdiff ! temporary scalar 1005 886 REAL(wp) :: zmax ! temporary scalar 1006 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt887 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt 1007 888 !!--------------------------------------------------------------------- 1008 889 ! 1009 ! IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 1010 ! 1011 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 890 ALLOCATE( zprt(jpi,jpj,jpk) ) 1012 891 ! 1013 892 IF(lwp) WRITE(numout,*) … … 1015 894 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 1016 895 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 896 897 ! compute position of the ice shelf grounding line 898 ! set bathy and isfdraft to 0 where grounded 899 IF ( ln_isfcav ) CALL zgr_isf_zspace 1017 900 1018 901 ! bathymetry in level (from bathy_meter) … … 1033 916 END DO 1034 917 918 ! Check compatibility between bathy and iceshelf draft 919 ! insure at least 2 wet level on the vertical under an ice shelf 920 ! compute misfdep and adjust isf draft if needed 921 IF ( ln_isfcav ) CALL zgr_isf_kspace 922 1035 923 ! Scale factors and depth at T- and W-points 1036 924 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1041 929 END DO 1042 930 1043 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf1044 IF ( ln_isfcav ) CALL zgr_isf1045 1046 931 ! Scale factors and depth at T- and W-points 1047 IF ( .NOT. ln_isfcav ) THEN1048 DO jj = 1, jpj1049 DO ji = 1, jpi1050 ik = mbathy(ji,jj)1051 IF( ik > 0 ) THEN ! ocean point only1052 ! max ocean level case1053 IF( ik == jpkm1 ) THEN1054 zdepwp = bathy(ji,jj)1055 ze3tp = bathy(ji,jj) - gdepw_1d(ik)1056 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )1057 e3t_0(ji,jj,ik ) = ze3tp1058 e3t_0(ji,jj,ik+1) = ze3tp1059 e3w_0(ji,jj,ik ) = ze3wp1060 e3w_0(ji,jj,ik+1) = ze3tp1061 gdepw_0(ji,jj,ik+1) = zdepwp1062 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp1063 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp1064 !1065 ELSE ! standard case1066 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)1067 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)1068 ENDIF1069 !gm Bug? check the gdepw_1d1070 ! ... on ik1071 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1072 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1073 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1074 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1075 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1076 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1077 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1078 ! ... on ik+11079 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1080 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1081 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1082 ENDIF1083 ENDIF1084 END DO1085 END DO1086 !1087 it = 01088 DO jj = 1, jpj1089 DO ji = 1, jpi1090 ik = mbathy(ji,jj)1091 IF( ik > 0 ) THEN ! ocean point only1092 e3tp (ji,jj) = e3t_0(ji,jj,ik)1093 e3wp (ji,jj) = e3w_0(ji,jj,ik)1094 ! test1095 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1096 IF( zdiff <= 0._wp .AND. lwp ) THEN1097 it = it + 11098 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1099 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1100 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1101 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1102 ENDIF1103 ENDIF1104 END DO1105 END DO1106 END IF1107 !1108 ! Scale factors and depth at U-, V-, UW and VW-points1109 DO jk = 1, jpk ! initialisation to z-scale factors1110 e3u_0 (:,:,jk) = e3t_1d(jk)1111 e3v_0 (:,:,jk) = e3t_1d(jk)1112 e3uw_0(:,:,jk) = e3w_1d(jk)1113 e3vw_0(:,:,jk) = e3w_1d(jk)1114 END DO1115 1116 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1117 DO jj = 1, jpjm11118 DO ji = 1, jpim1 ! vector opt.1119 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1120 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1121 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1122 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1123 END DO1124 END DO1125 END DO1126 IF ( ln_isfcav ) THEN1127 ! (ISF) define e3uw (adapted for 2 cells in the water column)1128 DO jj = 2, jpjm11129 DO ji = 2, jpim1 ! vector opt.1130 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1131 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1132 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1133 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1134 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1135 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1136 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1137 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1138 END DO1139 END DO1140 END IF1141 1142 CALL lbc_lnk('toto', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('toto', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1143 CALL lbc_lnk( 'toto',e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('toto', e3vw_0, 'V', 1._wp )1144 !1145 1146 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1147 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1148 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1149 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1150 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1151 END DO1152 1153 ! Scale factor at F-point1154 DO jk = 1, jpk ! initialisation to z-scale factors1155 e3f_0(:,:,jk) = e3t_1d(jk)1156 END DO1157 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1158 DO jj = 1, jpjm11159 DO ji = 1, jpim1 ! vector opt.1160 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1161 END DO1162 END DO1163 END DO1164 CALL lbc_lnk('toto', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1165 !1166 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1167 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1168 END DO1169 !!gm bug ? : must be a do loop with mj0,mj11170 !1171 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21172 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1173 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1174 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1175 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1176 1177 ! Control of the sign1178 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1179 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1180 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1181 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1182 1183 ! Compute gde3w_0 (vertical sum of e3w)1184 IF ( ln_isfcav ) THEN ! if cavity1185 WHERE( misfdep == 0 ) misfdep = 11186 DO jj = 1,jpj1187 DO ji = 1,jpi1188 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1189 DO jk = 2, misfdep(ji,jj)1190 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1191 END DO1192 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))1193 DO jk = misfdep(ji,jj) + 1, jpk1194 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1195 END DO1196 END DO1197 END DO1198 ELSE ! no cavity1199 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1200 DO jk = 2, jpk1201 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)1202 END DO1203 END IF1204 !1205 CALL wrk_dealloc( jpi,jpj,jpk, zprt )1206 !1207 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1208 !1209 END SUBROUTINE zgr_zps1210 1211 1212 SUBROUTINE zgr_isf1213 !!----------------------------------------------------------------------1214 !! *** ROUTINE zgr_isf ***1215 !!1216 !! ** Purpose : check the bathymetry in levels1217 !!1218 !! ** Method : THe water column have to contained at least 2 cells1219 !! Bathymetry and isfdraft are modified (dig/close) to respect1220 !! this criterion.1221 !!1222 !! ** Action : - test compatibility between isfdraft and bathy1223 !! - bathy and isfdraft are modified1224 !!----------------------------------------------------------------------1225 INTEGER :: ji, jj, jl, jk ! dummy loop indices1226 INTEGER :: ik, it ! temporary integers1227 INTEGER :: icompt, ibtest ! (ISF)1228 INTEGER :: ibtestim1, ibtestip1 ! (ISF)1229 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF)1230 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t1231 REAL(wp) :: zmax ! Maximum and minimum depth1232 REAL(wp) :: zbathydiff ! isf temporary scalar1233 REAL(wp) :: zrisfdepdiff ! isf temporary scalar1234 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1235 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t1236 REAL(wp) :: zdiff ! temporary scalar1237 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1238 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1239 !!---------------------------------------------------------------------1240 !1241 !! IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1242 !1243 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep)1244 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy )1245 1246 ! (ISF) compute misfdep1247 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 11248 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1249 END WHERE1250 1251 ! Compute misfdep for ocean points (i.e. first wet level)1252 ! find the first ocean level such that the first level thickness1253 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1254 ! e3t_0 is the reference level thickness1255 DO jk = 2, jpkm11256 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1257 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11258 END DO1259 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )1260 risfdep(:,:) = 0. ; misfdep(:,:) = 11261 END WHERE1262 1263 ! remove very shallow ice shelf (less than ~ 10m if 75L)1264 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1)1265 misfdep = 0; risfdep = 0.0_wp;1266 mbathy = 0; bathy = 0.0_wp;1267 END WHERE1268 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp)1269 misfdep = 0; risfdep = 0.0_wp;1270 mbathy = 0; bathy = 0.0_wp;1271 END WHERE1272 1273 ! 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 situation1274 icompt = 01275 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1276 DO jl = 1, 101277 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)1278 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin)1279 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1280 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1281 END WHERE1282 WHERE (mbathy(:,:) <= 0)1283 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1284 mbathy (:,:) = 0; bathy (:,:) = 0._wp1285 END WHERE1286 IF( lk_mpp ) THEN1287 zbathy(:,:) = FLOAT( misfdep(:,:) )1288 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1289 misfdep(:,:) = INT( zbathy(:,:) )1290 1291 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1292 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1293 1294 zbathy(:,:) = FLOAT( mbathy(:,:) )1295 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1296 mbathy(:,:) = INT( zbathy(:,:) )1297 ENDIF1298 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1299 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1300 misfdep(jpi,:) = misfdep( 2 ,:)1301 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1302 mbathy(jpi,:) = mbathy( 2 ,:)1303 ENDIF1304 1305 ! split last cell if possible (only where water column is 2 cell or less)1306 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).1307 IF ( .NOT. ln_iscpl) THEN1308 DO jk = jpkm1, 1, -11309 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1310 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1311 mbathy(:,:) = jk1312 bathy(:,:) = zmax1313 END WHERE1314 END DO1315 END IF1316 1317 ! split top cell if possible (only where water column is 2 cell or less)1318 DO jk = 2, jpkm11319 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1320 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1321 misfdep(:,:) = jk1322 risfdep(:,:) = zmax1323 END WHERE1324 END DO1325 1326 1327 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1328 DO jj = 1, jpj1329 DO ji = 1, jpi1330 ! find the minimum change option:1331 ! test bathy1332 IF (risfdep(ji,jj) > 1) THEN1333 IF ( .NOT. ln_iscpl ) THEN1334 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1335 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1336 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1337 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1338 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN1339 IF (zbathydiff <= zrisfdepdiff) THEN1340 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1341 mbathy(ji,jj)= mbathy(ji,jj) + 11342 ELSE1343 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1344 misfdep(ji,jj) = misfdep(ji,jj) - 11345 END IF1346 ENDIF1347 ELSE1348 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN1349 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1350 misfdep(ji,jj) = misfdep(ji,jj) - 11351 END IF1352 END IF1353 END IF1354 END DO1355 END DO1356 1357 ! At least 2 levels for water thickness at T, U, and V point.1358 DO jj = 1, jpj1359 DO ji = 1, jpi1360 ! find the minimum change option:1361 ! test bathy1362 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN1363 IF ( .NOT. ln_iscpl ) THEN1364 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &1365 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1366 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) &1367 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1368 IF (zbathydiff <= zrisfdepdiff) THEN1369 mbathy(ji,jj) = mbathy(ji,jj) + 11370 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1371 ELSE1372 misfdep(ji,jj)= misfdep(ji,jj) - 11373 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1374 END IF1375 ELSE1376 misfdep(ji,jj)= misfdep(ji,jj) - 11377 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1378 END IF1379 ENDIF1380 END DO1381 END DO1382 1383 ! point V mbathy(ji,jj) == misfdep(ji,jj+1)1384 DO jj = 1, jpjm11385 DO ji = 1, jpim11386 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN1387 IF ( .NOT. ln_iscpl ) THEN1388 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) &1389 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1390 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &1391 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1392 IF (zbathydiff <= zrisfdepdiff) THEN1393 mbathy(ji,jj) = mbathy(ji,jj) + 11394 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1395 ELSE1396 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11397 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1398 END IF1399 ELSE1400 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11401 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1402 END IF1403 ENDIF1404 END DO1405 END DO1406 1407 IF( lk_mpp ) THEN1408 zbathy(:,:) = FLOAT( misfdep(:,:) )1409 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1410 misfdep(:,:) = INT( zbathy(:,:) )1411 1412 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1413 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1414 1415 zbathy(:,:) = FLOAT( mbathy(:,:) )1416 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1417 mbathy(:,:) = INT( zbathy(:,:) )1418 ENDIF1419 ! point V misdep(ji,jj) == mbathy(ji,jj+1)1420 DO jj = 1, jpjm11421 DO ji = 1, jpim11422 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN1423 IF ( .NOT. ln_iscpl ) THEN1424 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &1425 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1426 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) &1427 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1428 IF (zbathydiff <= zrisfdepdiff) THEN1429 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11430 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1431 ELSE1432 misfdep(ji,jj) = misfdep(ji,jj) - 11433 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1434 END IF1435 ELSE1436 misfdep(ji,jj) = misfdep(ji,jj) - 11437 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1438 END IF1439 ENDIF1440 END DO1441 END DO1442 1443 1444 IF( lk_mpp ) THEN1445 zbathy(:,:) = FLOAT( misfdep(:,:) )1446 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1447 misfdep(:,:) = INT( zbathy(:,:) )1448 1449 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1450 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1451 1452 zbathy(:,:) = FLOAT( mbathy(:,:) )1453 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1454 mbathy(:,:) = INT( zbathy(:,:) )1455 ENDIF1456 1457 ! point U mbathy(ji,jj) == misfdep(ji,jj+1)1458 DO jj = 1, jpjm11459 DO ji = 1, jpim11460 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN1461 IF ( .NOT. ln_iscpl ) THEN1462 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &1463 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1464 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &1465 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1466 IF (zbathydiff <= zrisfdepdiff) THEN1467 mbathy(ji,jj) = mbathy(ji,jj) + 11468 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1469 ELSE1470 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11471 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1472 END IF1473 ELSE1474 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11475 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1476 ENDIF1477 ENDIF1478 ENDDO1479 ENDDO1480 1481 IF( lk_mpp ) THEN1482 zbathy(:,:) = FLOAT( misfdep(:,:) )1483 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1484 misfdep(:,:) = INT( zbathy(:,:) )1485 1486 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1487 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1488 1489 zbathy(:,:) = FLOAT( mbathy(:,:) )1490 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1491 mbathy(:,:) = INT( zbathy(:,:) )1492 ENDIF1493 1494 ! point U misfdep(ji,jj) == bathy(ji,jj+1)1495 DO jj = 1, jpjm11496 DO ji = 1, jpim11497 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN1498 IF ( .NOT. ln_iscpl ) THEN1499 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &1500 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1501 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) &1502 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1503 IF (zbathydiff <= zrisfdepdiff) THEN1504 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11505 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1506 ELSE1507 misfdep(ji,jj) = misfdep(ji ,jj) - 11508 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1509 END IF1510 ELSE1511 misfdep(ji,jj) = misfdep(ji ,jj) - 11512 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1513 ENDIF1514 ENDIF1515 ENDDO1516 ENDDO1517 1518 IF( lk_mpp ) THEN1519 zbathy(:,:) = FLOAT( misfdep(:,:) )1520 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1521 misfdep(:,:) = INT( zbathy(:,:) )1522 1523 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1524 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1525 1526 zbathy(:,:) = FLOAT( mbathy(:,:) )1527 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1528 mbathy(:,:) = INT( zbathy(:,:) )1529 ENDIF1530 END DO1531 ! end dig bathy/ice shelf to be compatible1532 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1533 DO jl = 1,201534 1535 ! remove single point "bay" on isf coast line in the ice shelf draft'1536 DO jk = 2, jpk1537 WHERE (misfdep==0) misfdep=jpk1538 zmask=0._wp1539 WHERE (misfdep <= jk) zmask=11540 DO jj = 2, jpjm11541 DO ji = 2, jpim11542 IF (misfdep(ji,jj) == jk) THEN1543 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1544 IF (ibtest <= 1) THEN1545 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11546 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk1547 END IF1548 END IF1549 END DO1550 END DO1551 END DO1552 WHERE (misfdep==jpk)1553 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp1554 END WHERE1555 IF( lk_mpp ) THEN1556 zbathy(:,:) = FLOAT( misfdep(:,:) )1557 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1558 misfdep(:,:) = INT( zbathy(:,:) )1559 1560 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1561 CALL lbc_lnk('toto', bathy, 'T', 1. )1562 1563 zbathy(:,:) = FLOAT( mbathy(:,:) )1564 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1565 mbathy(:,:) = INT( zbathy(:,:) )1566 ENDIF1567 1568 ! remove single point "bay" on bathy coast line beneath an ice shelf'1569 DO jk = jpk,1,-11570 zmask=0._wp1571 WHERE (mbathy >= jk ) zmask=11572 DO jj = 2, jpjm11573 DO ji = 2, jpim11574 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN1575 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1576 IF (ibtest <= 1) THEN1577 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11578 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 01579 END IF1580 END IF1581 END DO1582 END DO1583 END DO1584 WHERE (mbathy==0)1585 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp1586 END WHERE1587 IF( lk_mpp ) THEN1588 zbathy(:,:) = FLOAT( misfdep(:,:) )1589 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1590 misfdep(:,:) = INT( zbathy(:,:) )1591 1592 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1593 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1594 1595 zbathy(:,:) = FLOAT( mbathy(:,:) )1596 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1597 mbathy(:,:) = INT( zbathy(:,:) )1598 ENDIF1599 1600 ! fill hole in ice shelf1601 zmisfdep = misfdep1602 zrisfdep = risfdep1603 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk1604 DO jj = 2, jpjm11605 DO ji = 2, jpim11606 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1607 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1608 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk1609 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk1610 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk1611 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk1612 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1613 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN1614 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1615 END IF1616 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN1617 misfdep(ji,jj) = ibtest1618 risfdep(ji,jj) = gdepw_1d(ibtest)1619 ENDIF1620 ENDDO1621 ENDDO1622 1623 IF( lk_mpp ) THEN1624 zbathy(:,:) = FLOAT( misfdep(:,:) )1625 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1626 misfdep(:,:) = INT( zbathy(:,:) )1627 1628 CALL lbc_lnk( 'toto',risfdep, 'T', 1. )1629 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1630 1631 zbathy(:,:) = FLOAT( mbathy(:,:) )1632 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1633 mbathy(:,:) = INT( zbathy(:,:) )1634 ENDIF1635 !1636 !! fill hole in bathymetry1637 zmbathy (:,:)=mbathy (:,:)1638 DO jj = 2, jpjm11639 DO ji = 2, jpim11640 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1641 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1642 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 01643 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 01644 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 01645 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 01646 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1647 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN1648 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1649 END IF1650 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN1651 mbathy(ji,jj) = ibtest1652 bathy(ji,jj) = gdepw_1d(ibtest+1)1653 ENDIF1654 END DO1655 END DO1656 IF( lk_mpp ) THEN1657 zbathy(:,:) = FLOAT( misfdep(:,:) )1658 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1659 misfdep(:,:) = INT( zbathy(:,:) )1660 1661 CALL lbc_lnk( 'toto',risfdep, 'T', 1. )1662 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1663 1664 zbathy(:,:) = FLOAT( mbathy(:,:) )1665 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1666 mbathy(:,:) = INT( zbathy(:,:) )1667 ENDIF1668 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1669 DO jj = 1, jpjm11670 DO ji = 1, jpim11671 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN1672 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1673 END IF1674 END DO1675 END DO1676 IF( lk_mpp ) THEN1677 zbathy(:,:) = FLOAT( misfdep(:,:) )1678 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1679 misfdep(:,:) = INT( zbathy(:,:) )1680 1681 CALL lbc_lnk('toto', risfdep, 'T', 1. )1682 CALL lbc_lnk('toto', bathy, 'T', 1. )1683 1684 zbathy(:,:) = FLOAT( mbathy(:,:) )1685 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1686 mbathy(:,:) = INT( zbathy(:,:) )1687 ENDIF1688 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1689 DO jj = 1, jpjm11690 DO ji = 1, jpim11691 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN1692 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1693 END IF1694 END DO1695 END DO1696 IF( lk_mpp ) THEN1697 zbathy(:,:) = FLOAT( misfdep(:,:) )1698 CALL lbc_lnk('toto', zbathy, 'T', 1. )1699 misfdep(:,:) = INT( zbathy(:,:) )1700 1701 CALL lbc_lnk('toto', risfdep,'T', 1. )1702 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1703 1704 zbathy(:,:) = FLOAT( mbathy(:,:) )1705 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1706 mbathy(:,:) = INT( zbathy(:,:) )1707 ENDIF1708 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1709 DO jj = 1, jpjm11710 DO ji = 1, jpi1711 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN1712 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1713 END IF1714 END DO1715 END DO1716 IF( lk_mpp ) THEN1717 zbathy(:,:) = FLOAT( misfdep(:,:) )1718 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1719 misfdep(:,:) = INT( zbathy(:,:) )1720 1721 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1722 CALL lbc_lnk('toto', bathy, 'T', 1. )1723 1724 zbathy(:,:) = FLOAT( mbathy(:,:) )1725 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1726 mbathy(:,:) = INT( zbathy(:,:) )1727 ENDIF1728 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1729 DO jj = 1, jpjm11730 DO ji = 1, jpi1731 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN1732 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1733 END IF1734 END DO1735 END DO1736 IF( lk_mpp ) THEN1737 zbathy(:,:) = FLOAT( misfdep(:,:) )1738 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1739 misfdep(:,:) = INT( zbathy(:,:) )1740 1741 CALL lbc_lnk( 'toto',risfdep,'T', 1. )1742 CALL lbc_lnk( 'toto',bathy, 'T', 1. )1743 1744 zbathy(:,:) = FLOAT( mbathy(:,:) )1745 CALL lbc_lnk( 'toto',zbathy, 'T', 1. )1746 mbathy(:,:) = INT( zbathy(:,:) )1747 ENDIF1748 ! if not compatible after all check, mask T1749 DO jj = 1, jpj1750 DO ji = 1, jpi1751 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1752 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1753 END IF1754 END DO1755 END DO1756 1757 WHERE (mbathy(:,:) == 1)1758 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1759 END WHERE1760 END DO1761 ! end check compatibility ice shelf/bathy1762 ! remove very shallow ice shelf (less than ~ 10m if 75L)1763 WHERE (risfdep(:,:) <= 10._wp)1764 misfdep = 1; risfdep = 0.0_wp;1765 END WHERE1766 1767 IF( icompt == 0 ) THEN1768 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1769 ELSE1770 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1771 ENDIF1772 1773 ! compute scale factor and depth at T- and W- points1774 932 DO jj = 1, jpj 1775 933 DO ji = 1, jpi … … 1793 951 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1794 952 ENDIF 1795 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)1796 953 !gm Bug? check the gdepw_1d 1797 954 ! ... on ik … … 1799 956 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1800 957 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1801 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1802 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 958 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 959 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 960 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 961 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1803 962 ! ... on ik+1 1804 963 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1805 964 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 965 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1806 966 ENDIF 1807 967 ENDIF … … 1829 989 END DO 1830 990 ! 1831 ! (ISF) Definition of e3t, u, v, w for ISF case 1832 DO jj = 1, jpj 1833 DO ji = 1, jpi 1834 ik = misfdep(ji,jj) 1835 IF( ik > 1 ) THEN ! ice shelf point only 1836 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1837 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1838 !gm Bug? check the gdepw_0 1839 ! ... on ik 1840 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1841 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1842 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1843 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1844 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1845 1846 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1847 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1848 ENDIF 1849 ! ... on ik / ik-1 1850 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1851 gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik) 1852 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1853 e3w_0 (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2) 1854 gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1) 1855 ENDIF 1856 END DO 1857 END DO 1858 1859 it = 0 1860 DO jj = 1, jpj 1861 DO ji = 1, jpi 1862 ik = misfdep(ji,jj) 1863 IF( ik > 1 ) THEN ! ice shelf point only 1864 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1865 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1866 ! test 1867 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1868 IF( zdiff <= 0. .AND. lwp ) THEN 1869 it = it + 1 1870 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1871 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1872 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1873 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1874 ENDIF 1875 ENDIF 1876 END DO 1877 END DO 1878 1879 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1880 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1881 ! 1882 ! IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1883 ! 1884 END SUBROUTINE zgr_isf 1885 991 ! compute top scale factor if ice shelf 992 IF (ln_isfcav) CALL zps_isf 993 ! 994 ! Scale factors and depth at U-, V-, UW and VW-points 995 DO jk = 1, jpk ! initialisation to z-scale factors 996 e3u_0 (:,:,jk) = e3t_1d(jk) 997 e3v_0 (:,:,jk) = e3t_1d(jk) 998 e3uw_0(:,:,jk) = e3w_1d(jk) 999 e3vw_0(:,:,jk) = e3w_1d(jk) 1000 END DO 1001 1002 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1003 DO jj = 1, jpjm1 1004 DO ji = 1, jpim1 ! vector opt. 1005 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1006 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1007 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1008 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1009 END DO 1010 END DO 1011 END DO 1012 1013 ! update e3uw in case only 2 cells in the water column 1014 IF ( ln_isfcav ) CALL zps_isf_e3uv_w 1015 ! 1016 CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1017 CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 1018 ! 1019 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1020 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1021 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1022 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1023 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1024 END DO 1025 1026 ! Scale factor at F-point 1027 DO jk = 1, jpk ! initialisation to z-scale factors 1028 e3f_0(:,:,jk) = e3t_1d(jk) 1029 END DO 1030 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1031 DO jj = 1, jpjm1 1032 DO ji = 1, jpim1 ! vector opt. 1033 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1034 END DO 1035 END DO 1036 END DO 1037 CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1038 ! 1039 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1040 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1041 END DO 1042 !!gm bug ? : must be a do loop with mj0,mj1 1043 ! 1044 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1045 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1046 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1047 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1048 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1049 1050 ! Control of the sign 1051 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1052 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1053 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1054 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1055 ! 1056 ! if in the future gde3w_0 need to be compute, use the function defined in NEMO 1057 ! for now gde3w_0 computation is removed as not an output of domcfg 1058 1059 DEALLOCATE( zprt ) 1060 ! 1061 END SUBROUTINE zgr_zps 1886 1062 1887 1063 SUBROUTINE zgr_sco … … 1935 1111 REAL(wp) :: zrfact 1936 1112 ! 1937 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj21938 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat1113 REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1114 REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1939 1115 !! 1940 1116 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & … … 1942 1118 !!---------------------------------------------------------------------- 1943 1119 ! 1944 !! IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1945 ! 1946 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1120 ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) ) 1947 1121 ! 1948 1122 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 2024 1198 2025 1199 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 2026 CALL lbc_lnk( ' toto',zenv, 'T', 1._wp, 'no0' )1200 CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 2027 1201 ! 2028 1202 ! smooth the bathymetry (if required) … … 2088 1262 END DO 2089 1263 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 2090 CALL lbc_lnk( ' toto',zenv, 'T', 1._wp, 'no0' )1264 CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 2091 1265 ! ! ================ ! 2092 1266 END DO ! End loop ! … … 2132 1306 ! Apply lateral boundary condition 2133 1307 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 2134 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk(' toto', hbatu, 'U', 1._wp )1308 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp ) 2135 1309 DO jj = 1, jpj 2136 1310 DO ji = 1, jpi … … 2142 1316 END DO 2143 1317 END DO 2144 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk(' toto', hbatv, 'V', 1._wp )1318 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp ) 2145 1319 DO jj = 1, jpj 2146 1320 DO ji = 1, jpi … … 2151 1325 END DO 2152 1326 END DO 2153 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk(' toto', hbatf, 'F', 1._wp )1327 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp ) 2154 1328 DO jj = 1, jpj 2155 1329 DO ji = 1, jpi … … 2199 1373 ENDIF 2200 1374 2201 CALL lbc_lnk( ' toto',e3t_0 , 'T', 1._wp )2202 CALL lbc_lnk( ' toto',e3u_0 , 'U', 1._wp )2203 CALL lbc_lnk( ' toto',e3v_0 , 'V', 1._wp )2204 CALL lbc_lnk( ' toto',e3f_0 , 'F', 1._wp )2205 CALL lbc_lnk( ' toto',e3w_0 , 'W', 1._wp )2206 CALL lbc_lnk( ' toto',e3uw_0, 'U', 1._wp )2207 CALL lbc_lnk(' toto', e3vw_0, 'V', 1._wp )1375 CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp ) 1376 CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp ) 1377 CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp ) 1378 CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp ) 1379 CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp ) 1380 CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp ) 1381 CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 2208 1382 ! 2209 1383 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp … … 2214 1388 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2215 1389 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2216 2217 2218 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays)2219 !!gm and only that !!!!!2220 !!gm THIS should be removed from here !2221 gdept_n(:,:,:) = gdept_0(:,:,:)2222 gdepw_n(:,:,:) = gdepw_0(:,:,:)2223 gde3w_n(:,:,:) = gde3w_0(:,:,:)2224 e3t_n (:,:,:) = e3t_0 (:,:,:)2225 e3u_n (:,:,:) = e3u_0 (:,:,:)2226 e3v_n (:,:,:) = e3v_0 (:,:,:)2227 e3f_n (:,:,:) = e3f_0 (:,:,:)2228 e3w_n (:,:,:) = e3w_0 (:,:,:)2229 e3uw_n (:,:,:) = e3uw_0 (:,:,:)2230 e3vw_n (:,:,:) = e3vw_0 (:,:,:)2231 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine2232 !! gm end2233 1390 !! 2234 1391 ! HYBRID : … … 2236 1393 DO ji = 1, jpi 2237 1394 DO jk = 1, jpkm1 2238 IF( scobot(ji,jj) >= gdept_ n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )1395 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2239 1396 END DO 2240 1397 END DO … … 2246 1403 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2247 1404 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2248 & ' w ', MINVAL( gdepw_0(:,:,:) ) , '3w ' , MINVAL( gde3w_0(:,:,:) )1405 & ' w ', MINVAL( gdepw_0(:,:,:) ) 2249 1406 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2250 1407 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & … … 2253 1410 2254 1411 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2255 & ' w ', MAXVAL( gdepw_0(:,:,:) ) , '3w ' , MAXVAL( gde3w_0(:,:,:) )1412 & ' w ', MAXVAL( gdepw_0(:,:,:) ) 2256 1413 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2257 1414 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & … … 2298 1455 DO jk = 1, mbathy(ji,jj) 2299 1456 ! check coordinate is monotonically increasing 2300 IF (e3w_ n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN1457 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 2301 1458 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2302 1459 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2303 WRITE(numout,*) 'e3w',e3w_ n(ji,jj,:)2304 WRITE(numout,*) 'e3t',e3t_ n(ji,jj,:)1460 WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 1461 WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 2305 1462 CALL ctl_stop( ctmp1 ) 2306 1463 ENDIF 2307 1464 ! and check it has never gone negative 2308 IF( gdepw_ n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN1465 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 2309 1466 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2310 1467 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2311 WRITE(numout,*) 'gdepw',gdepw_ n(ji,jj,:)2312 WRITE(numout,*) 'gdept',gdept_ n(ji,jj,:)1468 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 1469 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 2313 1470 CALL ctl_stop( ctmp1 ) 2314 1471 ENDIF 2315 1472 ! and check it never exceeds the total depth 2316 IF( gdepw_ n(ji,jj,jk) > hbatt(ji,jj) ) THEN1473 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 2317 1474 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2318 1475 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2319 WRITE(numout,*) 'gdepw',gdepw_ n(ji,jj,:)1476 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 2320 1477 CALL ctl_stop( ctmp1 ) 2321 1478 ENDIF … … 2324 1481 DO jk = 1, mbathy(ji,jj)-1 2325 1482 ! and check it never exceeds the total depth 2326 IF( gdept_ n(ji,jj,jk) > hbatt(ji,jj) ) THEN1483 IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 2327 1484 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2328 1485 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2329 WRITE(numout,*) 'gdept',gdept_ n(ji,jj,:)1486 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 2330 1487 CALL ctl_stop( ctmp1 ) 2331 1488 ENDIF … … 2335 1492 END DO 2336 1493 ! 2337 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 2338 ! 2339 !!! IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1494 DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 2340 1495 ! 2341 1496 END SUBROUTINE zgr_sco … … 2358 1513 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2359 1514 ! 2360 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2361 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2362 !!---------------------------------------------------------------------- 2363 2364 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2365 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2366 2367 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1515 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3 1516 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1517 !!---------------------------------------------------------------------- 1518 1519 ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) ) 1520 ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) ) 1521 ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) ) 1522 1523 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp 2368 1524 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2369 1525 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp … … 2392 1548 z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 2393 1549 ! 2394 ! Coefficients for vertical depth as the sum of e3w scale factors2395 z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)2396 DO jk = 2, jpk2397 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)2398 END DO2399 !2400 1550 DO jk = 1, jpk 2401 1551 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) … … 2403 1553 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2404 1554 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2405 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2406 1555 END DO 2407 1556 ! … … 2448 1597 END DO 2449 1598 ! 2450 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3)2451 CALL wrk_dealloc( jpi,jpj,jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )1599 DEALLOCATE( z_gsigw3, z_gsigt3 ) 1600 DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2452 1601 ! 2453 1602 END SUBROUTINE s_sh94 … … 2476 1625 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2477 1626 ! 2478 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2479 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2480 !!---------------------------------------------------------------------- 2481 ! 2482 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2483 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2484 2485 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1627 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3 1628 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1629 !!---------------------------------------------------------------------- 1630 ! 1631 ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) ) 1632 ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk)) 1633 ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) ) 1634 1635 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp 2486 1636 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2487 1637 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp … … 2535 1685 z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 2536 1686 2537 ! Coefficients for vertical depth as the sum of e3w scale factors2538 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)2539 DO jk = 2, jpk2540 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)2541 END DO2542 2543 1687 DO jk = 1, jpk 2544 1688 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2545 1689 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2546 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2547 1690 END DO 2548 1691 … … 2608 1751 ENDDO 2609 1752 ! 2610 CALL lbc_lnk(' toto',e3t_0 ,'T',1.) ; CALL lbc_lnk('toto',e3u_0 ,'T',1.)2611 CALL lbc_lnk(' toto',e3v_0 ,'T',1.) ; CALL lbc_lnk('toto',e3f_0 ,'T',1.)2612 CALL lbc_lnk(' toto',e3w_0 ,'T',1.)2613 CALL lbc_lnk(' toto',e3uw_0,'T',1.) ; CALL lbc_lnk('toto',e3vw_0,'T',1.)2614 ! 2615 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3)2616 CALL wrk_dealloc( jpi,jpj,jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )1753 CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.) 1754 CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.) 1755 CALL lbc_lnk('domzgr',e3w_0 ,'T',1.) 1756 CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.) 1757 ! 1758 DEALLOCATE( z_gsigw3, z_gsigt3 ) 1759 DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2617 1760 ! 2618 1761 END SUBROUTINE s_sf12 … … 2631 1774 INTEGER :: ji, jj, jk ! dummy loop argument 2632 1775 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2633 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w2634 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw2635 !!---------------------------------------------------------------------- 2636 2637 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2638 CALL wrk_alloc( jpk, z_esigt, z_esigw)2639 2640 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp1776 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt 1777 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw 1778 !!---------------------------------------------------------------------- 1779 1780 ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) ) 1781 ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) 1782 1783 z_gsigw = 0._wp ; z_gsigt = 0._wp 2641 1784 z_esigt = 0._wp ; z_esigw = 0._wp 2642 1785 … … 2657 1800 z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 2658 1801 ! 2659 ! Coefficients for vertical depth as the sum of e3w scale factors2660 z_gsi3w(1) = 0.5_wp * z_esigw(1)2661 DO jk = 2, jpk2662 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)2663 END DO2664 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)2665 1802 DO jk = 1, jpk 2666 1803 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) … … 2668 1805 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2669 1806 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2670 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2671 END DO 2672 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1807 END DO 1808 2673 1809 DO jj = 1, jpj 2674 1810 DO ji = 1, jpi … … 2686 1822 END DO 2687 1823 ! 2688 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2689 CALL wrk_dealloc( jpk, z_esigt, z_esigw)1824 DEALLOCATE( z_gsigw, z_gsigt ) 1825 DEALLOCATE( z_esigt, z_esigw ) 2690 1826 ! 2691 1827 END SUBROUTINE s_tanh
Note: See TracChangeset
for help on using the changeset viewer.
