Changeset 6404 for branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2016-03-29T11:24:48+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6401 r6404 2 2 !!============================================================================== 3 3 !! *** MODULE domzgr *** 4 !! Ocean initialization : domain initialization4 !! Ocean domain : definition of the vertical coordinate system 5 5 !!============================================================================== 6 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate … … 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 20 21 !!---------------------------------------------------------------------- 21 22 … … 36 37 USE oce ! ocean variables 37 38 USE dom_oce ! ocean domain 39 USE wet_dry ! wetting and drying 38 40 USE closea ! closed seas 39 41 USE c1d ! 1D vertical configuration 42 ! 40 43 USE in_out_manager ! I/O manager 41 44 USE iom ! I/O library … … 73 76 74 77 !! * Substitutions 75 # include "domzgr_substitute.h90"76 78 # include "vectopt_loop_substitute.h90" 77 79 !!---------------------------------------------------------------------- … … 102 104 INTEGER :: ios 103 105 ! 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 105 107 !!---------------------------------------------------------------------- 106 108 ! … … 120 122 WRITE(numout,*) 'dom_zgr : vertical coordinate' 121 123 WRITE(numout,*) '~~~~~~~' 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 125 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 129 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 127 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 133 129 134 ioptio = 0 ! Check Vertical coordinate options … … 155 160 ! 156 161 IF( nprint == 1 .AND. lwp ) THEN 157 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )162 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 158 163 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 159 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )160 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &161 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &162 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &163 & ' w ', MINVAL(e3w_0(:,:,:) )164 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 165 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 166 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 167 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 168 & ' w ', MINVAL( e3w_0(:,:,:) ) 164 169 165 170 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 166 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )167 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &168 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &169 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&170 & ' w ', MAXVAL(e3w_0(:,:,:) )171 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 172 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 173 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 174 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 175 & ' w ', MAXVAL( e3w_0(:,:,:) ) 171 176 ENDIF 172 177 ! … … 384 389 !! - bathy : meter bathymetry (in meters) 385 390 !!---------------------------------------------------------------------- 386 INTEGER :: ji, jj, j l, jk ! dummy loop indices391 INTEGER :: ji, jj, jk ! dummy loop indices 387 392 INTEGER :: inum ! temporary logical unit 388 393 INTEGER :: ierror ! error flag … … 506 511 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 507 512 ! ! ===================== 508 IF( nn_cla == 0 ) THEN 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 514 END DO 513 ! 514 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 515 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 516 DO ji = mi0(ii0), mi1(ii1) 517 DO jj = mj0(ij0), mj1(ij1) 518 mbathy(ji,jj) = 15 515 519 END DO 516 IF(lwp) WRITE(numout,*)517 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0518 !519 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open520 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) = 12524 END DO520 END DO 521 IF(lwp) WRITE(numout,*) 522 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 523 ! 524 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 525 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 526 DO ji = mi0(ii0), mi1(ii1) 527 DO jj = mj0(ij0), mj1(ij1) 528 mbathy(ji,jj) = 12 525 529 END DO 526 IF(lwp) WRITE(numout,*)527 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0528 ENDIF530 END DO 531 IF(lwp) WRITE(numout,*) 532 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 529 533 ! 530 534 ENDIF … … 547 551 CALL iom_close( inum ) 548 552 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 553 554 ! set grounded point to 0 555 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 556 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 557 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 558 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 559 END WHERE 549 560 END IF 550 561 ! 551 562 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 552 ! 553 IF( nn_cla == 0 ) THEN 554 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 555 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 556 DO ji = mi0(ii0), mi1(ii1) 557 DO jj = mj0(ij0), mj1(ij1) 558 bathy(ji,jj) = 284._wp 559 END DO 563 ! 564 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 565 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 566 DO ji = mi0(ii0), mi1(ii1) 567 DO jj = mj0(ij0), mj1(ij1) 568 bathy(ji,jj) = 284._wp 560 569 END DO 561 IF(lwp) WRITE(numout,*)562 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0563 !564 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open565 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)566 DO ji = mi0(ii0), mi1(ii1)567 DO jj = mj0(ij0), mj1(ij1)568 bathy(ji,jj) = 137._wp569 END DO570 END DO 571 IF(lwp) WRITE(numout,*) 572 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 573 ! 574 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 575 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 576 DO ji = mi0(ii0), mi1(ii1) 577 DO jj = mj0(ij0), mj1(ij1) 578 bathy(ji,jj) = 137._wp 570 579 END DO 571 IF(lwp) WRITE(numout,*)572 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0573 ENDIF580 END DO 581 IF(lwp) WRITE(numout,*) 582 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 574 583 ! 575 584 ENDIF … … 586 595 ! 587 596 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 588 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin589 IF ( ln_isfcav ) THEN590 WHERE (bathy == risfdep)591 bathy = 0.0_wp ; risfdep = 0.0_wp592 END WHERE593 END IF594 ! end patch595 597 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 596 598 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 682 684 !! - update bathy : meter bathymetry (in meters) 683 685 !!---------------------------------------------------------------------- 684 !!685 686 INTEGER :: ji, jj, jl ! dummy loop indices 686 687 INTEGER :: icompt, ibtest, ikmax ! temporary integers 687 688 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 688 689 689 !!---------------------------------------------------------------------- 690 690 ! … … 783 783 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 784 784 ENDIF 785 786 IF( lwp .AND. nprint == 1 ) THEN ! control print787 WRITE(numout,*)788 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '789 WRITE(numout,*) ' ------------------'790 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )791 WRITE(numout,*)792 ENDIF793 785 ! 794 786 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 811 803 !! (min value = 1 over land) 812 804 !!---------------------------------------------------------------------- 813 !!814 805 INTEGER :: ji, jj ! dummy loop indices 815 806 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 843 834 END SUBROUTINE zgr_bot_level 844 835 845 SUBROUTINE zgr_top_level 846 !!---------------------------------------------------------------------- 847 !! *** ROUTINE zgr_bot_level *** 836 837 SUBROUTINE zgr_top_level 838 !!---------------------------------------------------------------------- 839 !! *** ROUTINE zgr_top_level *** 848 840 !! 849 841 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 855 847 !! (min value = 1) 856 848 !!---------------------------------------------------------------------- 857 !!858 849 INTEGER :: ji, jj ! dummy loop indices 859 850 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 889 880 END SUBROUTINE zgr_top_level 890 881 882 891 883 SUBROUTINE zgr_zco 892 884 !!---------------------------------------------------------------------- 893 885 !! *** ROUTINE zgr_zco *** 894 886 !! 895 !! ** Purpose : define the z-coordinate system887 !! ** Purpose : define the reference z-coordinate system 896 888 !! 897 889 !! ** Method : set 3D coord. arrays to reference 1D array … … 905 897 gdept_0(:,:,jk) = gdept_1d(jk) 906 898 gdepw_0(:,:,jk) = gdepw_1d(jk) 907 gde p3w_0(:,:,jk) = gdepw_1d(jk)899 gde3w_0(:,:,jk) = gdepw_1d(jk) 908 900 e3t_0(:,:,jk) = e3t_1d(jk) 909 901 e3u_0(:,:,jk) = e3t_1d(jk) … … 925 917 !! 926 918 !! ** Purpose : the depth and vertical scale factor in partial step 927 !! z-coordinate case919 !! reference z-coordinate case 928 920 !! 929 921 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 965 957 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 966 958 !!---------------------------------------------------------------------- 967 !!968 959 INTEGER :: ji, jj, jk ! dummy loop indices 969 960 INTEGER :: ik, it, ikb, ikt ! temporary integers 970 LOGICAL :: ll_print ! Allow control print for debugging971 961 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 972 962 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 973 REAL(wp) :: zmax ! Maximum depth974 963 REAL(wp) :: zdiff ! temporary scalar 975 REAL(wp) :: z refdep! temporary scalar964 REAL(wp) :: zmax ! temporary scalar 976 965 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 977 966 !!--------------------------------------------------------------------- … … 979 968 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 980 969 ! 981 CALL wrk_alloc( jpi, jpj, jpk,zprt )970 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 982 971 ! 983 972 IF(lwp) WRITE(numout,*) … … 985 974 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 986 975 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 987 988 ll_print = .FALSE. ! Local variable for debugging989 990 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth991 WRITE(numout,*)992 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'993 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )994 ENDIF995 996 976 997 977 ! bathymetry in level (from bathy_meter) … … 1012 992 END DO 1013 993 1014 IF ( ln_isfcav ) CALL zgr_isf1015 1016 994 ! Scale factors and depth at T- and W-points 1017 995 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1021 999 e3w_0 (:,:,jk) = e3w_1d (jk) 1022 1000 END DO 1001 1002 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1003 IF ( ln_isfcav ) CALL zgr_isf 1004 1005 ! Scale factors and depth at T- and W-points 1006 IF ( .NOT. ln_isfcav ) THEN 1007 DO jj = 1, jpj 1008 DO ji = 1, jpi 1009 ik = mbathy(ji,jj) 1010 IF( ik > 0 ) THEN ! ocean point only 1011 ! max ocean level case 1012 IF( ik == jpkm1 ) THEN 1013 zdepwp = bathy(ji,jj) 1014 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1015 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1016 e3t_0(ji,jj,ik ) = ze3tp 1017 e3t_0(ji,jj,ik+1) = ze3tp 1018 e3w_0(ji,jj,ik ) = ze3wp 1019 e3w_0(ji,jj,ik+1) = ze3tp 1020 gdepw_0(ji,jj,ik+1) = zdepwp 1021 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1022 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1023 ! 1024 ELSE ! standard case 1025 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1026 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1027 ENDIF 1028 !gm Bug? check the gdepw_1d 1029 ! ... on ik 1030 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1031 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1032 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1033 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1034 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1035 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1036 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1037 ! ... on ik+1 1038 e3w_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 1039 e3t_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 1040 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1041 ENDIF 1042 ENDIF 1043 END DO 1044 END DO 1045 ! 1046 it = 0 1047 DO jj = 1, jpj 1048 DO ji = 1, jpi 1049 ik = mbathy(ji,jj) 1050 IF( ik > 0 ) THEN ! ocean point only 1051 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1052 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1053 ! test 1054 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1055 IF( zdiff <= 0._wp .AND. lwp ) THEN 1056 it = it + 1 1057 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1058 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1059 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1060 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1061 ENDIF 1062 ENDIF 1063 END DO 1064 END DO 1065 END IF 1066 ! 1067 ! Scale factors and depth at U-, V-, UW and VW-points 1068 DO jk = 1, jpk ! initialisation to z-scale factors 1069 e3u_0 (:,:,jk) = e3t_1d(jk) 1070 e3v_0 (:,:,jk) = e3t_1d(jk) 1071 e3uw_0(:,:,jk) = e3w_1d(jk) 1072 e3vw_0(:,:,jk) = e3w_1d(jk) 1073 END DO 1074 1075 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1076 DO jj = 1, jpjm1 1077 DO ji = 1, fs_jpim1 ! vector opt. 1078 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1079 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1080 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1081 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1082 END DO 1083 END DO 1084 END DO 1085 IF ( ln_isfcav ) THEN 1086 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1087 DO jj = 2, jpjm1 1088 DO ji = 2, fs_jpim1 ! vector opt. 1089 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1090 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1091 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1092 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1093 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1094 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1095 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1096 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1097 END DO 1098 END DO 1099 END IF 1100 1101 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1102 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1103 ! 1104 1105 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1106 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1107 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1108 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1109 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1110 END DO 1111 1112 ! Scale factor at F-point 1113 DO jk = 1, jpk ! initialisation to z-scale factors 1114 e3f_0(:,:,jk) = e3t_1d(jk) 1115 END DO 1116 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1117 DO jj = 1, jpjm1 1118 DO ji = 1, fs_jpim1 ! vector opt. 1119 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1120 END DO 1121 END DO 1122 END DO 1123 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1124 ! 1125 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1126 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1127 END DO 1128 !!gm bug ? : must be a do loop with mj0,mj1 1023 1129 ! 1130 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1131 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1132 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1133 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1134 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1135 1136 ! Control of the sign 1137 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1138 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1139 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1140 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1141 1142 ! Compute gde3w_0 (vertical sum of e3w) 1143 IF ( ln_isfcav ) THEN ! if cavity 1144 WHERE( misfdep == 0 ) misfdep = 1 1145 DO jj = 1,jpj 1146 DO ji = 1,jpi 1147 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1148 DO jk = 2, misfdep(ji,jj) 1149 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1150 END DO 1151 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)) 1152 DO jk = misfdep(ji,jj) + 1, jpk 1153 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1154 END DO 1155 END DO 1156 END DO 1157 ELSE ! no cavity 1158 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1159 DO jk = 2, jpk 1160 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1161 END DO 1162 END IF 1163 ! 1164 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1165 ! 1166 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1167 ! 1168 END SUBROUTINE zgr_zps 1169 1170 1171 SUBROUTINE zgr_isf 1172 !!---------------------------------------------------------------------- 1173 !! *** ROUTINE zgr_isf *** 1174 !! 1175 !! ** Purpose : check the bathymetry in levels 1176 !! 1177 !! ** Method : THe water column have to contained at least 2 cells 1178 !! Bathymetry and isfdraft are modified (dig/close) to respect 1179 !! this criterion. 1180 !! 1181 !! ** Action : - test compatibility between isfdraft and bathy 1182 !! - bathy and isfdraft are modified 1183 !!---------------------------------------------------------------------- 1184 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1185 INTEGER :: ik, it ! temporary integers 1186 INTEGER :: icompt, ibtest ! (ISF) 1187 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1188 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1189 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1190 REAL(wp) :: zmax ! Maximum and minimum depth 1191 REAL(wp) :: zbathydiff ! isf temporary scalar 1192 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1193 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1194 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1195 REAL(wp) :: zdiff ! temporary scalar 1196 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1197 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1198 !!--------------------------------------------------------------------- 1199 ! 1200 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1201 ! 1202 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1203 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1204 1205 ! (ISF) compute misfdep 1206 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1207 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1208 END WHERE 1209 1210 ! Compute misfdep for ocean points (i.e. first wet level) 1211 ! find the first ocean level such that the first level thickness 1212 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1213 ! e3t_0 is the reference level thickness 1214 DO jk = 2, jpkm1 1215 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1216 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1217 END DO 1218 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1219 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1220 END WHERE 1221 1222 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1223 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1224 misfdep = 0; risfdep = 0.0_wp; 1225 mbathy = 0; bathy = 0.0_wp; 1226 END WHERE 1227 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1228 misfdep = 0; risfdep = 0.0_wp; 1229 mbathy = 0; bathy = 0.0_wp; 1230 END WHERE 1231 1232 ! 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 1233 icompt = 0 1234 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1235 DO jl = 1, 10 1236 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1237 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1238 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1239 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1240 END WHERE 1241 WHERE (mbathy(:,:) <= 0) 1242 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1243 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1244 END WHERE 1245 IF( lk_mpp ) THEN 1246 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1247 CALL lbc_lnk( zbathy, 'T', 1. ) 1248 misfdep(:,:) = INT( zbathy(:,:) ) 1249 1250 CALL lbc_lnk( risfdep,'T', 1. ) 1251 CALL lbc_lnk( bathy, 'T', 1. ) 1252 1253 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1254 CALL lbc_lnk( zbathy, 'T', 1. ) 1255 mbathy(:,:) = INT( zbathy(:,:) ) 1256 ENDIF 1257 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1258 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1259 misfdep(jpi,:) = misfdep( 2 ,:) 1260 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1261 mbathy(jpi,:) = mbathy( 2 ,:) 1262 ENDIF 1263 1264 ! split last cell if possible (only where water column is 2 cell or less) 1265 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1266 IF ( .NOT. ln_iscpl) THEN 1267 DO jk = jpkm1, 1, -1 1268 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1269 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1270 mbathy(:,:) = jk 1271 bathy(:,:) = zmax 1272 END WHERE 1273 END DO 1274 END IF 1275 1276 ! split top cell if possible (only where water column is 2 cell or less) 1277 DO jk = 2, jpkm1 1278 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1279 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1280 misfdep(:,:) = jk 1281 risfdep(:,:) = zmax 1282 END WHERE 1283 END DO 1284 1285 1286 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1287 DO jj = 1, jpj 1288 DO ji = 1, jpi 1289 ! find the minimum change option: 1290 ! test bathy 1291 IF (risfdep(ji,jj) > 1) THEN 1292 IF ( .NOT. ln_iscpl ) THEN 1293 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1294 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1295 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1296 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1297 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1298 IF (zbathydiff <= zrisfdepdiff) THEN 1299 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1300 mbathy(ji,jj)= mbathy(ji,jj) + 1 1301 ELSE 1302 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1303 misfdep(ji,jj) = misfdep(ji,jj) - 1 1304 END IF 1305 ENDIF 1306 ELSE 1307 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1308 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1309 misfdep(ji,jj) = misfdep(ji,jj) - 1 1310 END IF 1311 END IF 1312 END IF 1313 END DO 1314 END DO 1315 1316 ! At least 2 levels for water thickness at T, U, and V point. 1317 DO jj = 1, jpj 1318 DO ji = 1, jpi 1319 ! find the minimum change option: 1320 ! test bathy 1321 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1322 IF ( .NOT. ln_iscpl ) THEN 1323 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1324 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1325 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1326 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1327 IF (zbathydiff <= zrisfdepdiff) THEN 1328 mbathy(ji,jj) = mbathy(ji,jj) + 1 1329 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1330 ELSE 1331 misfdep(ji,jj)= misfdep(ji,jj) - 1 1332 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1333 END IF 1334 ELSE 1335 misfdep(ji,jj)= misfdep(ji,jj) - 1 1336 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1337 END IF 1338 ENDIF 1339 END DO 1340 END DO 1341 1342 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1343 DO jj = 1, jpjm1 1344 DO ji = 1, jpim1 1345 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1346 IF ( .NOT. ln_iscpl ) THEN 1347 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1348 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1349 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1350 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1351 IF (zbathydiff <= zrisfdepdiff) THEN 1352 mbathy(ji,jj) = mbathy(ji,jj) + 1 1353 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1354 ELSE 1355 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1356 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1357 END IF 1358 ELSE 1359 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1360 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1361 END IF 1362 ENDIF 1363 END DO 1364 END DO 1365 1366 IF( lk_mpp ) THEN 1367 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1368 CALL lbc_lnk( zbathy, 'T', 1. ) 1369 misfdep(:,:) = INT( zbathy(:,:) ) 1370 1371 CALL lbc_lnk( risfdep,'T', 1. ) 1372 CALL lbc_lnk( bathy, 'T', 1. ) 1373 1374 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1375 CALL lbc_lnk( zbathy, 'T', 1. ) 1376 mbathy(:,:) = INT( zbathy(:,:) ) 1377 ENDIF 1378 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1379 DO jj = 1, jpjm1 1380 DO ji = 1, jpim1 1381 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1382 IF ( .NOT. ln_iscpl ) THEN 1383 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1384 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1385 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1386 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1387 IF (zbathydiff <= zrisfdepdiff) THEN 1388 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1389 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1390 ELSE 1391 misfdep(ji,jj) = misfdep(ji,jj) - 1 1392 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1393 END IF 1394 ELSE 1395 misfdep(ji,jj) = misfdep(ji,jj) - 1 1396 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1397 END IF 1398 ENDIF 1399 END DO 1400 END DO 1401 1402 1403 IF( lk_mpp ) THEN 1404 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1405 CALL lbc_lnk( zbathy, 'T', 1. ) 1406 misfdep(:,:) = INT( zbathy(:,:) ) 1407 1408 CALL lbc_lnk( risfdep,'T', 1. ) 1409 CALL lbc_lnk( bathy, 'T', 1. ) 1410 1411 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1412 CALL lbc_lnk( zbathy, 'T', 1. ) 1413 mbathy(:,:) = INT( zbathy(:,:) ) 1414 ENDIF 1415 1416 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1417 DO jj = 1, jpjm1 1418 DO ji = 1, jpim1 1419 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1420 IF ( .NOT. ln_iscpl ) THEN 1421 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1422 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1423 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1424 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1425 IF (zbathydiff <= zrisfdepdiff) THEN 1426 mbathy(ji,jj) = mbathy(ji,jj) + 1 1427 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1428 ELSE 1429 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1430 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1431 END IF 1432 ELSE 1433 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1434 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1435 ENDIF 1436 ENDIF 1437 ENDDO 1438 ENDDO 1439 1440 IF( lk_mpp ) THEN 1441 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1442 CALL lbc_lnk( zbathy, 'T', 1. ) 1443 misfdep(:,:) = INT( zbathy(:,:) ) 1444 1445 CALL lbc_lnk( risfdep,'T', 1. ) 1446 CALL lbc_lnk( bathy, 'T', 1. ) 1447 1448 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1449 CALL lbc_lnk( zbathy, 'T', 1. ) 1450 mbathy(:,:) = INT( zbathy(:,:) ) 1451 ENDIF 1452 1453 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1454 DO jj = 1, jpjm1 1455 DO ji = 1, jpim1 1456 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1457 IF ( .NOT. ln_iscpl ) THEN 1458 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1459 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1460 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1461 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1462 IF (zbathydiff <= zrisfdepdiff) THEN 1463 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1464 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1465 ELSE 1466 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1467 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1468 END IF 1469 ELSE 1470 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1471 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1472 ENDIF 1473 ENDIF 1474 ENDDO 1475 ENDDO 1476 1477 IF( lk_mpp ) THEN 1478 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1479 CALL lbc_lnk( zbathy, 'T', 1. ) 1480 misfdep(:,:) = INT( zbathy(:,:) ) 1481 1482 CALL lbc_lnk( risfdep,'T', 1. ) 1483 CALL lbc_lnk( bathy, 'T', 1. ) 1484 1485 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1486 CALL lbc_lnk( zbathy, 'T', 1. ) 1487 mbathy(:,:) = INT( zbathy(:,:) ) 1488 ENDIF 1489 END DO 1490 ! end dig bathy/ice shelf to be compatible 1491 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1492 DO jl = 1,20 1493 1494 ! remove single point "bay" on isf coast line in the ice shelf draft' 1495 DO jk = 2, jpk 1496 WHERE (misfdep==0) misfdep=jpk 1497 zmask=0._wp 1498 WHERE (misfdep <= jk) zmask=1 1499 DO jj = 2, jpjm1 1500 DO ji = 2, jpim1 1501 IF (misfdep(ji,jj) == jk) THEN 1502 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1503 IF (ibtest <= 1) THEN 1504 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1505 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1506 END IF 1507 END IF 1508 END DO 1509 END DO 1510 END DO 1511 WHERE (misfdep==jpk) 1512 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1513 END WHERE 1514 IF( lk_mpp ) THEN 1515 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1516 CALL lbc_lnk( zbathy, 'T', 1. ) 1517 misfdep(:,:) = INT( zbathy(:,:) ) 1518 1519 CALL lbc_lnk( risfdep,'T', 1. ) 1520 CALL lbc_lnk( bathy, 'T', 1. ) 1521 1522 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1523 CALL lbc_lnk( zbathy, 'T', 1. ) 1524 mbathy(:,:) = INT( zbathy(:,:) ) 1525 ENDIF 1526 1527 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1528 DO jk = jpk,1,-1 1529 zmask=0._wp 1530 WHERE (mbathy >= jk ) zmask=1 1531 DO jj = 2, jpjm1 1532 DO ji = 2, jpim1 1533 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1534 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1535 IF (ibtest <= 1) THEN 1536 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1537 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1538 END IF 1539 END IF 1540 END DO 1541 END DO 1542 END DO 1543 WHERE (mbathy==0) 1544 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1545 END WHERE 1546 IF( lk_mpp ) THEN 1547 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1548 CALL lbc_lnk( zbathy, 'T', 1. ) 1549 misfdep(:,:) = INT( zbathy(:,:) ) 1550 1551 CALL lbc_lnk( risfdep,'T', 1. ) 1552 CALL lbc_lnk( bathy, 'T', 1. ) 1553 1554 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1555 CALL lbc_lnk( zbathy, 'T', 1. ) 1556 mbathy(:,:) = INT( zbathy(:,:) ) 1557 ENDIF 1558 1559 ! fill hole in ice shelf 1560 zmisfdep = misfdep 1561 zrisfdep = risfdep 1562 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1563 DO jj = 2, jpjm1 1564 DO ji = 2, jpim1 1565 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1566 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1567 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1568 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1569 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1570 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1571 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1572 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1573 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1574 END IF 1575 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1576 misfdep(ji,jj) = ibtest 1577 risfdep(ji,jj) = gdepw_1d(ibtest) 1578 ENDIF 1579 ENDDO 1580 ENDDO 1581 1582 IF( lk_mpp ) THEN 1583 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1584 CALL lbc_lnk( zbathy, 'T', 1. ) 1585 misfdep(:,:) = INT( zbathy(:,:) ) 1586 1587 CALL lbc_lnk( risfdep, 'T', 1. ) 1588 CALL lbc_lnk( bathy, 'T', 1. ) 1589 1590 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1591 CALL lbc_lnk( zbathy, 'T', 1. ) 1592 mbathy(:,:) = INT( zbathy(:,:) ) 1593 ENDIF 1594 ! 1595 !! fill hole in bathymetry 1596 zmbathy (:,:)=mbathy (:,:) 1597 DO jj = 2, jpjm1 1598 DO ji = 2, jpim1 1599 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1600 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1601 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1602 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1603 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1604 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1605 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1606 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1607 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1608 END IF 1609 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1610 mbathy(ji,jj) = ibtest 1611 bathy(ji,jj) = gdepw_1d(ibtest+1) 1612 ENDIF 1613 END DO 1614 END DO 1615 IF( lk_mpp ) THEN 1616 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1617 CALL lbc_lnk( zbathy, 'T', 1. ) 1618 misfdep(:,:) = INT( zbathy(:,:) ) 1619 1620 CALL lbc_lnk( risfdep, 'T', 1. ) 1621 CALL lbc_lnk( bathy, 'T', 1. ) 1622 1623 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1624 CALL lbc_lnk( zbathy, 'T', 1. ) 1625 mbathy(:,:) = INT( zbathy(:,:) ) 1626 ENDIF 1627 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1628 DO jj = 1, jpjm1 1629 DO ji = 1, jpim1 1630 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1631 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1632 END IF 1633 END DO 1634 END DO 1635 IF( lk_mpp ) THEN 1636 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1637 CALL lbc_lnk( zbathy, 'T', 1. ) 1638 misfdep(:,:) = INT( zbathy(:,:) ) 1639 1640 CALL lbc_lnk( risfdep, 'T', 1. ) 1641 CALL lbc_lnk( bathy, 'T', 1. ) 1642 1643 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1644 CALL lbc_lnk( zbathy, 'T', 1. ) 1645 mbathy(:,:) = INT( zbathy(:,:) ) 1646 ENDIF 1647 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1648 DO jj = 1, jpjm1 1649 DO ji = 1, jpim1 1650 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1651 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1652 END IF 1653 END DO 1654 END DO 1655 IF( lk_mpp ) THEN 1656 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1657 CALL lbc_lnk( zbathy, 'T', 1. ) 1658 misfdep(:,:) = INT( zbathy(:,:) ) 1659 1660 CALL lbc_lnk( risfdep,'T', 1. ) 1661 CALL lbc_lnk( bathy, 'T', 1. ) 1662 1663 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1664 CALL lbc_lnk( zbathy, 'T', 1. ) 1665 mbathy(:,:) = INT( zbathy(:,:) ) 1666 ENDIF 1667 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1668 DO jj = 1, jpjm1 1669 DO ji = 1, jpi 1670 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1671 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1672 END IF 1673 END DO 1674 END DO 1675 IF( lk_mpp ) THEN 1676 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1677 CALL lbc_lnk( zbathy, 'T', 1. ) 1678 misfdep(:,:) = INT( zbathy(:,:) ) 1679 1680 CALL lbc_lnk( risfdep,'T', 1. ) 1681 CALL lbc_lnk( bathy, 'T', 1. ) 1682 1683 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1684 CALL lbc_lnk( zbathy, 'T', 1. ) 1685 mbathy(:,:) = INT( zbathy(:,:) ) 1686 ENDIF 1687 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1688 DO jj = 1, jpjm1 1689 DO ji = 1, jpi 1690 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1691 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1692 END IF 1693 END DO 1694 END DO 1695 IF( lk_mpp ) THEN 1696 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1697 CALL lbc_lnk( zbathy, 'T', 1. ) 1698 misfdep(:,:) = INT( zbathy(:,:) ) 1699 1700 CALL lbc_lnk( risfdep,'T', 1. ) 1701 CALL lbc_lnk( bathy, 'T', 1. ) 1702 1703 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1704 CALL lbc_lnk( zbathy, 'T', 1. ) 1705 mbathy(:,:) = INT( zbathy(:,:) ) 1706 ENDIF 1707 ! if not compatible after all check, mask T 1708 DO jj = 1, jpj 1709 DO ji = 1, jpi 1710 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1711 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1712 END IF 1713 END DO 1714 END DO 1715 1716 WHERE (mbathy(:,:) == 1) 1717 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1718 END WHERE 1719 END DO 1720 ! end check compatibility ice shelf/bathy 1721 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1722 WHERE (risfdep(:,:) <= 10._wp) 1723 misfdep = 1; risfdep = 0.0_wp; 1724 END WHERE 1725 1726 IF( icompt == 0 ) THEN 1727 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1728 ELSE 1729 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1730 ENDIF 1731 1732 ! compute scale factor and depth at T- and W- points 1024 1733 DO jj = 1, jpj 1025 1734 DO ji = 1, jpi … … 1043 1752 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1044 1753 ENDIF 1754 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1045 1755 !gm Bug? check the gdepw_1d 1046 1756 ! ... on ik … … 1048 1758 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1049 1759 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1050 e3t_0(ji,jj,ik) = e3t_1d(ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1051 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1052 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1053 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1760 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1761 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1054 1762 ! ... on ik+1 1055 e3w_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 1056 e3t_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 1057 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1763 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1764 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1058 1765 ENDIF 1059 1766 ENDIF … … 1081 1788 END DO 1082 1789 ! 1083 IF ( ln_isfcav ) THEN1084 1790 ! (ISF) Definition of e3t, u, v, w for ISF case 1085 1086 1087 1088 1089 1090 1791 DO jj = 1, jpj 1792 DO ji = 1, jpi 1793 ik = misfdep(ji,jj) 1794 IF( ik > 1 ) THEN ! ice shelf point only 1795 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1796 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1091 1797 !gm Bug? check the gdepw_0 1092 1093 1094 1095 1096 e3t_0(ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)1097 e3w_0(ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)1098 1099 1100 e3w_0(ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)1101 1102 1103 e3w_0 (ji,jj,ik ) =2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1104 1798 ! ... on ik 1799 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1800 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1801 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1802 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1803 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1804 1805 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1806 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1807 ENDIF 1808 ! ... on ik / ik-1 1809 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1810 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1105 1811 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1106 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1812 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1813 ENDIF 1814 END DO 1815 END DO 1816 1817 it = 0 1818 DO jj = 1, jpj 1819 DO ji = 1, jpi 1820 ik = misfdep(ji,jj) 1821 IF( ik > 1 ) THEN ! ice shelf point only 1822 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1823 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1824 ! test 1825 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1826 IF( zdiff <= 0. .AND. lwp ) THEN 1827 it = it + 1 1828 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1829 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1830 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1831 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1107 1832 ENDIF 1108 END DO1833 ENDIF 1109 1834 END DO 1110 !1111 it = 01112 DO jj = 1, jpj1113 DO ji = 1, jpi1114 ik = misfdep(ji,jj)1115 IF( ik > 1 ) THEN ! ice shelf point only1116 e3tp (ji,jj) = e3t_0(ji,jj,ik )1117 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1118 ! test1119 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1120 IF( zdiff <= 0. .AND. lwp ) THEN1121 it = it + 11122 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1123 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1124 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1125 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1126 ENDIF1127 ENDIF1128 END DO1129 END DO1130 END IF1131 ! END (ISF)1132 1133 ! Scale factors and depth at U-, V-, UW and VW-points1134 DO jk = 1, jpk ! initialisation to z-scale factors1135 e3u_0 (:,:,jk) = e3t_1d(jk)1136 e3v_0 (:,:,jk) = e3t_1d(jk)1137 e3uw_0(:,:,jk) = e3w_1d(jk)1138 e3vw_0(:,:,jk) = e3w_1d(jk)1139 END DO1140 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1141 DO jj = 1, jpjm11142 DO ji = 1, fs_jpim1 ! vector opt.1143 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1144 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1145 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1146 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1147 END DO1148 END DO1149 END DO1150 IF ( ln_isfcav ) THEN1151 ! (ISF) define e3uw (adapted for 2 cells in the water column)1152 DO jj = 2, jpjm11153 DO ji = 2, fs_jpim1 ! vector opt.1154 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1155 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1156 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1157 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1158 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1159 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1160 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1161 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1162 END DO1163 END DO1164 END IF1165 1166 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1167 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1168 !1169 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1170 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1171 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1172 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1173 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1174 END DO1175 1176 ! Scale factor at F-point1177 DO jk = 1, jpk ! initialisation to z-scale factors1178 e3f_0(:,:,jk) = e3t_1d(jk)1179 END DO1180 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1181 DO jj = 1, jpjm11182 DO ji = 1, fs_jpim1 ! vector opt.1183 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1184 END DO1185 END DO1186 END DO1187 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1188 !1189 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1190 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1191 END DO1192 !!gm bug ? : must be a do loop with mj0,mj11193 !1194 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21195 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1196 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1197 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1198 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1199 1200 ! Control of the sign1201 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1202 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1203 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1204 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1205 1206 ! Compute gdep3w_0 (vertical sum of e3w)1207 IF ( ln_isfcav ) THEN ! if cavity1208 WHERE (misfdep == 0) misfdep = 11209 DO jj = 1,jpj1210 DO ji = 1,jpi1211 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1212 DO jk = 2, misfdep(ji,jj)1213 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1214 END DO1215 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1216 DO jk = misfdep(ji,jj) + 1, jpk1217 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1218 END DO1219 END DO1220 END DO1221 ELSE ! no cavity1222 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1223 DO jk = 2, jpk1224 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1225 END DO1226 END IF1227 ! ! ================= !1228 IF(lwp .AND. ll_print) THEN ! Control print !1229 ! ! ================= !1230 DO jj = 1,jpj1231 DO ji = 1, jpi1232 ik = MAX( mbathy(ji,jj), 1 )1233 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1234 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1235 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1236 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1237 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1238 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1239 END DO1240 END DO1241 WRITE(numout,*)1242 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1243 WRITE(numout,*)1244 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1245 WRITE(numout,*)1246 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1247 WRITE(numout,*)1248 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1249 WRITE(numout,*)1250 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1251 WRITE(numout,*)1252 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1253 ENDIF1254 !1255 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1256 !1257 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1258 !1259 END SUBROUTINE zgr_zps1260 1261 SUBROUTINE zgr_isf1262 !!----------------------------------------------------------------------1263 !! *** ROUTINE zgr_isf ***1264 !!1265 !! ** Purpose : check the bathymetry in levels1266 !!1267 !! ** Method : THe water column have to contained at least 2 cells1268 !! Bathymetry and isfdraft are modified (dig/close) to respect1269 !! this criterion.1270 !!1271 !!1272 !! ** Action : - test compatibility between isfdraft and bathy1273 !! - bathy and isfdraft are modified1274 !!----------------------------------------------------------------------1275 !!1276 INTEGER :: ji, jj, jk, jl ! dummy loop indices1277 INTEGER :: ik, it ! temporary integers1278 INTEGER :: id, jd, nprocd1279 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1280 LOGICAL :: ll_print ! Allow control print for debugging1281 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1282 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1283 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1284 REAL(wp) :: zdiff ! temporary scalar1285 REAL(wp) :: zrefdep ! temporary scalar1286 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1287 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1288 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1289 !!---------------------------------------------------------------------1290 !1291 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1292 !1293 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)1294 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )1295 1296 1297 ! (ISF) compute misfdep1298 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11299 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1300 END WHERE1301 1302 ! Compute misfdep for ocean points (i.e. first wet level)1303 ! find the first ocean level such that the first level thickness1304 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1305 ! e3t_0 is the reference level thickness1306 DO jk = 2, jpkm11307 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1308 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11309 1835 END DO 1310 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1311 risfdep(:,:) = 0. ; misfdep(:,:) = 11312 END WHERE1313 1314 ! 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 situation1315 icompt = 01316 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1317 DO jl = 1, 101318 WHERE (bathy(:,:) .EQ. risfdep(:,:) )1319 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1320 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1321 END WHERE1322 WHERE (mbathy(:,:) <= 0)1323 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1324 mbathy (:,:) = 0; bathy (:,:) = 0._wp1325 ENDWHERE1326 IF( lk_mpp ) THEN1327 zbathy(:,:) = FLOAT( misfdep(:,:) )1328 CALL lbc_lnk( zbathy, 'T', 1. )1329 misfdep(:,:) = INT( zbathy(:,:) )1330 CALL lbc_lnk( risfdep, 'T', 1. )1331 CALL lbc_lnk( bathy, 'T', 1. )1332 zbathy(:,:) = FLOAT( mbathy(:,:) )1333 CALL lbc_lnk( zbathy, 'T', 1. )1334 mbathy(:,:) = INT( zbathy(:,:) )1335 ENDIF1336 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1337 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1338 misfdep(jpi,:) = misfdep( 2 ,:)1339 ENDIF1340 1341 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1342 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1343 mbathy(jpi,:) = mbathy( 2 ,:)1344 ENDIF1345 1346 ! split last cell if possible (only where water column is 2 cell or less)1347 DO jk = jpkm1, 1, -11348 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1349 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1350 mbathy(:,:) = jk1351 bathy(:,:) = zmax1352 END WHERE1353 END DO1354 1355 ! split top cell if possible (only where water column is 2 cell or less)1356 DO jk = 2, jpkm11357 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1358 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1359 misfdep(:,:) = jk1360 risfdep(:,:) = zmax1361 END WHERE1362 END DO1363 1364 1365 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1366 DO jj = 1, jpj1367 DO ji = 1, jpi1368 ! find the minimum change option:1369 ! test bathy1370 IF (risfdep(ji,jj) .GT. 1) THEN1371 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1372 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1373 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1374 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1375 1376 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN1377 IF (zbathydiff .LE. zrisfdepdiff) THEN1378 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1379 mbathy(ji,jj)= mbathy(ji,jj) + 11380 ELSE1381 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1382 misfdep(ji,jj) = misfdep(ji,jj) - 11383 END IF1384 END IF1385 END IF1386 END DO1387 END DO1388 1389 ! At least 2 levels for water thickness at T, U, and V point.1390 DO jj = 1, jpj1391 DO ji = 1, jpi1392 ! find the minimum change option:1393 ! test bathy1394 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1395 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1396 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1397 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1398 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1399 IF (zbathydiff .LE. zrisfdepdiff) THEN1400 mbathy(ji,jj) = mbathy(ji,jj) + 11401 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1402 ELSE1403 misfdep(ji,jj)= misfdep(ji,jj) - 11404 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1405 END IF1406 ENDIF1407 END DO1408 END DO1409 1410 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1411 DO jj = 1, jpjm11412 DO ji = 1, jpim11413 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1414 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1415 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1416 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1417 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1418 IF (zbathydiff .LE. zrisfdepdiff) THEN1419 mbathy(ji,jj) = mbathy(ji,jj) + 11420 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1421 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1422 ELSE1423 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11424 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1425 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1426 END IF1427 ENDIF1428 END DO1429 END DO1430 1431 IF( lk_mpp ) THEN1432 zbathy(:,:) = FLOAT( misfdep(:,:) )1433 CALL lbc_lnk( zbathy, 'T', 1. )1434 misfdep(:,:) = INT( zbathy(:,:) )1435 CALL lbc_lnk( risfdep, 'T', 1. )1436 CALL lbc_lnk( bathy, 'T', 1. )1437 zbathy(:,:) = FLOAT( mbathy(:,:) )1438 CALL lbc_lnk( zbathy, 'T', 1. )1439 mbathy(:,:) = INT( zbathy(:,:) )1440 ENDIF1441 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1442 DO jj = 1, jpjm11443 DO ji = 1, jpim11444 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1445 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1446 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1447 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1448 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1449 IF (zbathydiff .LE. zrisfdepdiff) THEN1450 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11451 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1452 ELSE1453 misfdep(ji,jj) = misfdep(ji,jj) - 11454 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1455 END IF1456 ENDIF1457 END DO1458 END DO1459 1460 1461 IF( lk_mpp ) THEN1462 zbathy(:,:) = FLOAT( misfdep(:,:) )1463 CALL lbc_lnk( zbathy, 'T', 1. )1464 misfdep(:,:) = INT( zbathy(:,:) )1465 CALL lbc_lnk( risfdep, 'T', 1. )1466 CALL lbc_lnk( bathy, 'T', 1. )1467 zbathy(:,:) = FLOAT( mbathy(:,:) )1468 CALL lbc_lnk( zbathy, 'T', 1. )1469 mbathy(:,:) = INT( zbathy(:,:) )1470 ENDIF1471 1472 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1473 DO jj = 1, jpjm11474 DO ji = 1, jpim11475 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1476 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1477 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1478 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1479 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1480 IF (zbathydiff .LE. zrisfdepdiff) THEN1481 mbathy(ji,jj) = mbathy(ji,jj) + 11482 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1483 ELSE1484 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11485 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1486 END IF1487 ENDIF1488 ENDDO1489 ENDDO1490 1491 IF( lk_mpp ) THEN1492 zbathy(:,:) = FLOAT( misfdep(:,:) )1493 CALL lbc_lnk( zbathy, 'T', 1. )1494 misfdep(:,:) = INT( zbathy(:,:) )1495 CALL lbc_lnk( risfdep, 'T', 1. )1496 CALL lbc_lnk( bathy, 'T', 1. )1497 zbathy(:,:) = FLOAT( mbathy(:,:) )1498 CALL lbc_lnk( zbathy, 'T', 1. )1499 mbathy(:,:) = INT( zbathy(:,:) )1500 ENDIF1501 1502 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1503 DO jj = 1, jpjm11504 DO ji = 1, jpim11505 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1506 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1507 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1508 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1509 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1510 IF (zbathydiff .LE. zrisfdepdiff) THEN1511 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11512 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1513 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1514 ELSE1515 misfdep(ji,jj) = misfdep(ji ,jj) - 11516 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1517 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1518 END IF1519 ENDIF1520 ENDDO1521 ENDDO1522 1523 IF( lk_mpp ) THEN1524 zbathy(:,:) = FLOAT( misfdep(:,:) )1525 CALL lbc_lnk( zbathy, 'T', 1. )1526 misfdep(:,:) = INT( zbathy(:,:) )1527 CALL lbc_lnk( risfdep, 'T', 1. )1528 CALL lbc_lnk( bathy, 'T', 1. )1529 zbathy(:,:) = FLOAT( mbathy(:,:) )1530 CALL lbc_lnk( zbathy, 'T', 1. )1531 mbathy(:,:) = INT( zbathy(:,:) )1532 ENDIF1533 END DO1534 ! end dig bathy/ice shelf to be compatible1535 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1536 DO jl = 1,201537 1538 ! remove single point "bay" on isf coast line in the ice shelf draft'1539 DO jk = 2, jpk1540 WHERE (misfdep==0) misfdep=jpk1541 zmask=01542 WHERE (misfdep .LE. jk) zmask=11543 DO jj = 2, jpjm11544 DO ji = 2, jpim11545 IF (misfdep(ji,jj) .EQ. jk) THEN1546 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1547 IF (ibtest .LE. 1) THEN1548 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11549 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1550 END IF1551 END IF1552 END DO1553 END DO1554 END DO1555 WHERE (misfdep==jpk)1556 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1557 END WHERE1558 IF( lk_mpp ) THEN1559 zbathy(:,:) = FLOAT( misfdep(:,:) )1560 CALL lbc_lnk( zbathy, 'T', 1. )1561 misfdep(:,:) = INT( zbathy(:,:) )1562 CALL lbc_lnk( risfdep, 'T', 1. )1563 CALL lbc_lnk( bathy, 'T', 1. )1564 zbathy(:,:) = FLOAT( mbathy(:,:) )1565 CALL lbc_lnk( zbathy, 'T', 1. )1566 mbathy(:,:) = INT( zbathy(:,:) )1567 ENDIF1568 1569 ! remove single point "bay" on bathy coast line beneath an ice shelf'1570 DO jk = jpk,1,-11571 zmask=01572 WHERE (mbathy .GE. jk ) zmask=11573 DO jj = 2, jpjm11574 DO ji = 2, jpim11575 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1576 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1577 IF (ibtest .LE. 1) THEN1578 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11579 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01580 END IF1581 END IF1582 END DO1583 END DO1584 END DO1585 WHERE (mbathy==0)1586 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1587 END WHERE1588 IF( lk_mpp ) THEN1589 zbathy(:,:) = FLOAT( misfdep(:,:) )1590 CALL lbc_lnk( zbathy, 'T', 1. )1591 misfdep(:,:) = INT( zbathy(:,:) )1592 CALL lbc_lnk( risfdep, 'T', 1. )1593 CALL lbc_lnk( bathy, 'T', 1. )1594 zbathy(:,:) = FLOAT( mbathy(:,:) )1595 CALL lbc_lnk( zbathy, 'T', 1. )1596 mbathy(:,:) = INT( zbathy(:,:) )1597 ENDIF1598 1599 ! fill hole in ice shelf1600 zmisfdep = misfdep1601 zrisfdep = risfdep1602 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1603 DO jj = 2, jpjm11604 DO ji = 2, jpim11605 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1606 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1607 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1608 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1609 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1610 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1611 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1612 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1613 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1614 END IF1615 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1616 misfdep(ji,jj) = ibtest1617 risfdep(ji,jj) = gdepw_1d(ibtest)1618 ENDIF1619 ENDDO1620 ENDDO1621 1622 IF( lk_mpp ) THEN1623 zbathy(:,:) = FLOAT( misfdep(:,:) )1624 CALL lbc_lnk( zbathy, 'T', 1. )1625 misfdep(:,:) = INT( zbathy(:,:) )1626 CALL lbc_lnk( risfdep, 'T', 1. )1627 CALL lbc_lnk( bathy, 'T', 1. )1628 zbathy(:,:) = FLOAT( mbathy(:,:) )1629 CALL lbc_lnk( zbathy, 'T', 1. )1630 mbathy(:,:) = INT( zbathy(:,:) )1631 ENDIF1632 !1633 !! fill hole in bathymetry1634 zmbathy (:,:)=mbathy (:,:)1635 DO jj = 2, jpjm11636 DO ji = 2, jpim11637 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1638 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1639 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1640 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01641 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01642 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01643 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1644 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1645 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1646 END IF1647 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1648 mbathy(ji,jj) = ibtest1649 bathy(ji,jj) = gdepw_1d(ibtest+1)1650 ENDIF1651 END DO1652 END DO1653 IF( lk_mpp ) THEN1654 zbathy(:,:) = FLOAT( misfdep(:,:) )1655 CALL lbc_lnk( zbathy, 'T', 1. )1656 misfdep(:,:) = INT( zbathy(:,:) )1657 CALL lbc_lnk( risfdep, 'T', 1. )1658 CALL lbc_lnk( bathy, 'T', 1. )1659 zbathy(:,:) = FLOAT( mbathy(:,:) )1660 CALL lbc_lnk( zbathy, 'T', 1. )1661 mbathy(:,:) = INT( zbathy(:,:) )1662 ENDIF1663 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1664 DO jj = 1, jpjm11665 DO ji = 1, jpim11666 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1667 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1668 END IF1669 END DO1670 END DO1671 IF( lk_mpp ) THEN1672 zbathy(:,:) = FLOAT( misfdep(:,:) )1673 CALL lbc_lnk( zbathy, 'T', 1. )1674 misfdep(:,:) = INT( zbathy(:,:) )1675 CALL lbc_lnk( risfdep, 'T', 1. )1676 CALL lbc_lnk( bathy, 'T', 1. )1677 zbathy(:,:) = FLOAT( mbathy(:,:) )1678 CALL lbc_lnk( zbathy, 'T', 1. )1679 mbathy(:,:) = INT( zbathy(:,:) )1680 ENDIF1681 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1682 DO jj = 1, jpjm11683 DO ji = 1, jpim11684 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1685 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1686 END IF1687 END DO1688 END DO1689 IF( lk_mpp ) THEN1690 zbathy(:,:) = FLOAT( misfdep(:,:) )1691 CALL lbc_lnk( zbathy, 'T', 1. )1692 misfdep(:,:) = INT( zbathy(:,:) )1693 CALL lbc_lnk( risfdep, 'T', 1. )1694 CALL lbc_lnk( bathy, 'T', 1. )1695 zbathy(:,:) = FLOAT( mbathy(:,:) )1696 CALL lbc_lnk( zbathy, 'T', 1. )1697 mbathy(:,:) = INT( zbathy(:,:) )1698 ENDIF1699 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1700 DO jj = 1, jpjm11701 DO ji = 1, jpi1702 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1703 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1704 END IF1705 END DO1706 END DO1707 IF( lk_mpp ) THEN1708 zbathy(:,:) = FLOAT( misfdep(:,:) )1709 CALL lbc_lnk( zbathy, 'T', 1. )1710 misfdep(:,:) = INT( zbathy(:,:) )1711 CALL lbc_lnk( risfdep, 'T', 1. )1712 CALL lbc_lnk( bathy, 'T', 1. )1713 zbathy(:,:) = FLOAT( mbathy(:,:) )1714 CALL lbc_lnk( zbathy, 'T', 1. )1715 mbathy(:,:) = INT( zbathy(:,:) )1716 ENDIF1717 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1718 DO jj = 1, jpjm11719 DO ji = 1, jpi1720 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1721 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1722 END IF1723 END DO1724 END DO1725 IF( lk_mpp ) THEN1726 zbathy(:,:) = FLOAT( misfdep(:,:) )1727 CALL lbc_lnk( zbathy, 'T', 1. )1728 misfdep(:,:) = INT( zbathy(:,:) )1729 CALL lbc_lnk( risfdep, 'T', 1. )1730 CALL lbc_lnk( bathy, 'T', 1. )1731 zbathy(:,:) = FLOAT( mbathy(:,:) )1732 CALL lbc_lnk( zbathy, 'T', 1. )1733 mbathy(:,:) = INT( zbathy(:,:) )1734 ENDIF1735 ! if not compatible after all check, mask T1736 DO jj = 1, jpj1737 DO ji = 1, jpi1738 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1739 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1740 END IF1741 END DO1742 END DO1743 1744 WHERE (mbathy(:,:) == 1)1745 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1746 END WHERE1747 END DO1748 ! end check compatibility ice shelf/bathy1749 ! remove very shallow ice shelf (less than ~ 10m if 75L)1750 WHERE (misfdep(:,:) <= 5)1751 misfdep = 1; risfdep = 0.0_wp;1752 END WHERE1753 1754 IF( icompt == 0 ) THEN1755 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1756 ELSE1757 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1758 ENDIF1759 1836 1760 1837 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1761 1838 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1762 1763 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1764 1765 END SUBROUTINE 1839 ! 1840 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1841 ! 1842 END SUBROUTINE zgr_isf 1843 1766 1844 1767 1845 SUBROUTINE zgr_sco … … 1809 1887 !! 1810 1888 !!---------------------------------------------------------------------- 1811 !1812 1889 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1813 1890 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1818 1895 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1819 1896 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1820 1897 !! 1821 1898 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1822 1899 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1823 1900 !!---------------------------------------------------------------------- 1824 1901 ! 1825 1902 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1826 1903 ! 1827 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1904 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1828 1905 ! 1829 1906 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1870 1947 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1871 1948 1872 DO jj = 1, jpj 1873 DO ji = 1, jpi 1874 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1875 END DO 1876 END DO 1949 IF( .NOT.ln_wd ) THEN 1950 DO jj = 1, jpj 1951 DO ji = 1, jpi 1952 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1953 END DO 1954 END DO 1955 END IF 1877 1956 ! ! ============================= 1878 1957 ! ! Define the envelop bathymetry (hbatt) … … 1881 1960 zenv(:,:) = bathy(:,:) 1882 1961 ! 1962 IF( .NOT.ln_wd ) THEN 1883 1963 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1884 DO jj = 1, jpj 1885 DO ji = 1, jpi 1886 IF( bathy(ji,jj) == 0._wp ) THEN 1887 iip1 = MIN( ji+1, jpi ) 1888 ijp1 = MIN( jj+1, jpj ) 1889 iim1 = MAX( ji-1, 1 ) 1890 ijm1 = MAX( jj-1, 1 ) 1891 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1892 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1893 zenv(ji,jj) = rn_sbot_min 1894 ENDIF 1895 ENDIF 1896 END DO 1897 END DO 1964 DO jj = 1, jpj 1965 DO ji = 1, jpi 1966 IF( bathy(ji,jj) == 0._wp ) THEN 1967 iip1 = MIN( ji+1, jpi ) 1968 ijp1 = MIN( jj+1, jpj ) 1969 iim1 = MAX( ji-1, 1 ) 1970 ijm1 = MAX( jj-1, 1 ) 1971 !!gm BUG fix see ticket #1617 1972 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1973 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1974 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) & 1975 & zenv(ji,jj) = rn_sbot_min 1976 !!gm 1977 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1978 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1979 !!gm zenv(ji,jj) = rn_sbot_min 1980 !!gm ENDIF 1981 !!gm end 1982 ENDIF 1983 END DO 1984 END DO 1985 END IF 1986 1898 1987 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1899 1988 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1983 2072 ENDIF 1984 2073 ! 1985 IF(lwp) THEN ! Control print1986 WRITE(numout,*)1987 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1988 WRITE(numout,*)1989 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1990 IF( nprint == 1 ) THEN1991 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1992 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1993 ENDIF1994 ENDIF1995 1996 2074 ! ! ============================== 1997 2075 ! ! hbatu, hbatv, hbatf fields … … 1999 2077 IF(lwp) THEN 2000 2078 WRITE(numout,*) 2001 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2079 IF( .NOT.ln_wd ) THEN 2080 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2081 ELSE 2082 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2083 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2084 ENDIF 2002 2085 ENDIF 2003 2086 hbatu(:,:) = rn_sbot_min … … 2012 2095 END DO 2013 2096 END DO 2097 2098 IF( ln_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points 2099 DO jj = 1, jpj 2100 DO ji = 1, jpi 2101 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2102 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2103 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2104 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2105 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2106 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2107 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2108 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2109 END DO 2110 END DO 2111 END IF 2014 2112 ! 2015 2113 ! Apply lateral boundary condition … … 2019 2117 DO ji = 1, jpi 2020 2118 IF( hbatu(ji,jj) == 0._wp ) THEN 2119 !No worries about the following line when ln_wd == .true. 2021 2120 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2022 2121 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 2044 2143 2045 2144 !!bug: key_helsinki a verifer 2046 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2047 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2048 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2049 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2145 IF( .NOT.ln_wd ) THEN 2146 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2147 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2148 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2149 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2150 END IF 2050 2151 2051 2152 IF( nprint == 1 .AND. lwp ) THEN … … 2088 2189 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2089 2190 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2090 2091 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2092 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2093 ! 2094 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2095 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2096 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2097 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2098 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2099 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2100 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2191 ! 2192 IF( .NOT.ln_wd ) THEN 2193 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2194 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2195 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2196 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2197 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2198 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2199 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2200 END IF 2101 2201 2102 2202 #if defined key_agrif 2103 ! Ensure meaningful vertical scale factors in ghost lines/columns 2104 IF( .NOT. Agrif_Root() ) THEN 2105 ! 2106 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2107 e3u_0(1,:,:) = e3u_0(2,:,:) 2108 ENDIF 2109 ! 2110 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2111 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2112 ENDIF 2113 ! 2114 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2115 e3v_0(:,1,:) = e3v_0(:,2,:) 2116 ENDIF 2117 ! 2118 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2119 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2120 ENDIF 2121 ! 2122 ENDIF 2203 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2204 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2205 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2206 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2207 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2208 ENDIF 2123 2209 #endif 2124 2210 2125 fsdept(:,:,:) = gdept_0 (:,:,:) 2126 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2127 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2128 fse3t(:,:,:) = e3t_0(:,:,:) 2129 fse3u(:,:,:) = e3u_0(:,:,:) 2130 fse3v(:,:,:) = e3v_0(:,:,:) 2131 fse3f(:,:,:) = e3f_0(:,:,:) 2132 fse3w(:,:,:) = e3w_0(:,:,:) 2133 fse3uw(:,:,:) = e3uw_0(:,:,:) 2134 fse3vw(:,:,:) = e3vw_0(:,:,:) 2211 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2212 !!gm and only that !!!!! 2213 !!gm THIS should be removed from here ! 2214 gdept_n(:,:,:) = gdept_0(:,:,:) 2215 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2216 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2217 e3t_n (:,:,:) = e3t_0 (:,:,:) 2218 e3u_n (:,:,:) = e3u_0 (:,:,:) 2219 e3v_n (:,:,:) = e3v_0 (:,:,:) 2220 e3f_n (:,:,:) = e3f_0 (:,:,:) 2221 e3w_n (:,:,:) = e3w_0 (:,:,:) 2222 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2223 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2224 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2225 !! gm end 2135 2226 !! 2136 2227 ! HYBRID : … … 2138 2229 DO ji = 1, jpi 2139 2230 DO jk = 1, jpkm1 2140 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2141 END DO 2142 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2231 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2232 END DO 2233 IF( ln_wd ) THEN 2234 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2235 mbathy(ji,jj) = 0 2236 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2237 mbathy(ji,jj) = 2 2238 ENDIF 2239 ELSE 2240 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2241 ENDIF 2143 2242 END DO 2144 2243 END DO … … 2149 2248 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2150 2249 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2151 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2152 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2153 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2154 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2250 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2251 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2252 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2253 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2155 2254 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2156 2255 2157 2256 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2158 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2159 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2160 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2161 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2257 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2258 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2259 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2260 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2162 2261 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2163 2262 ENDIF … … 2191 2290 END DO 2192 2291 ENDIF 2193 2194 !================================================================================2195 ! check the coordinate makes sense2196 !================================================================================2292 ! 2293 !================================================================================ 2294 ! check the coordinate makes sense 2295 !================================================================================ 2197 2296 DO ji = 1, jpi 2198 2297 DO jj = 1, jpj 2199 2298 ! 2200 2299 IF( hbatt(ji,jj) > 0._wp) THEN 2201 2300 DO jk = 1, mbathy(ji,jj) 2202 2301 ! check coordinate is monotonically increasing 2203 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2302 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2204 2303 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2205 2304 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2206 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2207 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2305 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2306 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2208 2307 CALL ctl_stop( ctmp1 ) 2209 2308 ENDIF 2210 2309 ! and check it has never gone negative 2211 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2310 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2212 2311 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2213 2312 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2214 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2215 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2313 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2314 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2216 2315 CALL ctl_stop( ctmp1 ) 2217 2316 ENDIF 2218 2317 ! and check it never exceeds the total depth 2219 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2318 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2220 2319 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2221 2320 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2222 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2321 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2223 2322 CALL ctl_stop( ctmp1 ) 2224 2323 ENDIF 2225 2324 END DO 2226 2325 ! 2227 2326 DO jk = 1, mbathy(ji,jj)-1 2228 2327 ! and check it never exceeds the total depth 2229 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2328 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2230 2329 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2231 2330 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2232 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2331 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2233 2332 CALL ctl_stop( ctmp1 ) 2234 2333 ENDIF 2235 2334 END DO 2236 2237 2335 ENDIF 2238 2239 2336 END DO 2240 2337 END DO … … 2246 2343 END SUBROUTINE zgr_sco 2247 2344 2248 !!====================================================================== 2345 2249 2346 SUBROUTINE s_sh94() 2250 2251 2347 !!---------------------------------------------------------------------- 2252 2348 !! *** ROUTINE s_sh94 *** … … 2259 2355 !! Reference : Song and Haidvogel 1994. 2260 2356 !!---------------------------------------------------------------------- 2261 !2262 2357 INTEGER :: ji, jj, jk ! dummy loop argument 2263 2358 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2359 REAL(wp) :: ztmpu, ztmpv, ztmpf 2360 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2264 2361 ! 2265 2362 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2266 2363 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2267 2268 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2269 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2364 !!---------------------------------------------------------------------- 2365 2366 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2367 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2270 2368 2271 2369 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2273 2371 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2274 2372 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2275 2373 ! 2276 2374 DO ji = 1, jpi 2277 2375 DO jj = 1, jpj 2278 2376 ! 2279 2377 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2280 2378 DO jk = 1, jpk … … 2305 2403 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2306 2404 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2307 gdept_0 2308 gdepw_0 2309 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2405 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2406 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2407 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2310 2408 END DO 2311 2409 ! … … 2315 2413 DO ji = 1, jpim1 2316 2414 DO jj = 1, jpjm1 2415 ! extended for Wetting/Drying case 2416 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2417 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2418 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2419 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2420 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2421 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2422 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2317 2423 DO jk = 1, jpk 2318 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2319 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2320 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2321 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2322 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2323 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2324 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2325 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2326 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2327 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2328 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2424 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2425 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2426 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2427 ELSE 2428 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2429 & / ztmpu 2430 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2431 & / ztmpu 2432 END IF 2433 2434 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2435 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2436 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2437 ELSE 2438 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2439 & / ztmpv 2440 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2441 & / ztmpv 2442 END IF 2443 2444 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2445 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj ,jk) + z_esigt3(ji+1,jj ,jk) & 2446 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2447 ELSE 2448 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2449 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2450 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2451 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2452 END IF 2453 2329 2454 ! 2330 2455 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 2339 2464 END DO 2340 2465 END DO 2341 2342 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2343 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2344 2466 ! 2467 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2468 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2469 ! 2345 2470 END SUBROUTINE s_sh94 2346 2471 2472 2347 2473 SUBROUTINE s_sf12 2348 2349 2474 !!---------------------------------------------------------------------- 2350 2475 !! *** ROUTINE s_sf12 *** … … 2362 2487 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2363 2488 !!---------------------------------------------------------------------- 2364 !2365 2489 INTEGER :: ji, jj, jk ! dummy loop argument 2366 2490 REAL(wp) :: zsmth ! smoothing around critical depth 2367 2491 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2492 REAL(wp) :: ztmpu, ztmpv, ztmpf 2493 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2368 2494 ! 2369 2495 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2370 2496 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2371 2497 !!---------------------------------------------------------------------- 2372 2498 ! 2373 2499 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2433 2559 2434 2560 DO jk = 1, jpk 2435 gdept_0 2436 gdepw_0 2437 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2561 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2562 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2563 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2438 2564 END DO 2439 2565 … … 2444 2570 DO jj=1,jpj-1 2445 2571 2446 DO jk = 1, jpk 2447 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 2448 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2449 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 2450 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2451 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 2452 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 2453 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2454 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 2455 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2456 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 2457 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2572 ! extend to suit for Wetting/Drying case 2573 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2574 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2575 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2576 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2577 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2578 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2579 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2580 DO jk = 1, jpk 2581 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2582 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2583 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2584 ELSE 2585 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2586 & / ztmpu 2587 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2588 & / ztmpu 2589 END IF 2590 2591 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2592 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2593 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2594 ELSE 2595 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2596 & / ztmpv 2597 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2598 & / ztmpv 2599 END IF 2600 2601 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2602 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2603 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2604 ELSE 2605 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2606 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2607 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2608 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2609 END IF 2610 2611 ! Code prior to wetting and drying option (for reference) 2612 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2613 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2614 ! 2615 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2616 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2617 ! 2618 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2619 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2620 ! 2621 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2622 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2623 ! 2624 !z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2625 ! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2626 ! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2627 ! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2628 ! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2458 2629 2459 2630 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) … … 2462 2633 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2463 2634 ! 2464 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2635 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2465 2636 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2466 2637 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2475 2646 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2476 2647 ! 2477 ! ! ============= 2478 2479 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2480 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2481 2648 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2649 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2650 ! 2482 2651 END SUBROUTINE s_sf12 2483 2652 2653 2484 2654 SUBROUTINE s_tanh() 2485 2486 2655 !!---------------------------------------------------------------------- 2487 2656 !! *** ROUTINE s_tanh*** … … 2493 2662 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2494 2663 !!---------------------------------------------------------------------- 2495 2496 INTEGER :: ji, jj, jk ! dummy loop argument 2664 INTEGER :: ji, jj, jk ! dummy loop argument 2497 2665 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2498 2499 2666 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2500 2667 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2501 2502 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2503 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2668 !!---------------------------------------------------------------------- 2669 2670 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2671 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2504 2672 2505 2673 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2531 2699 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2532 2700 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2533 gdept_0 2534 gdepw_0 2535 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2701 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2702 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2703 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2536 2704 END DO 2537 2705 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2550 2718 END DO 2551 2719 END DO 2552 2553 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2554 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2555 2720 ! 2721 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2722 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2723 ! 2556 2724 END SUBROUTINE s_tanh 2725 2557 2726 2558 2727 FUNCTION fssig( pk ) RESULT( pf ) … … 2626 2795 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2627 2796 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2628 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2629 REAL(wp), INTENT(in ) :: pzs ! surface box depth2630 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2631 REAL(wp) :: za1,za2,za3 ! local variables2632 REAL(wp) :: zn1,zn2 ! local variables2633 REAL(wp) :: za,zb,zx ! local variables2634 integer :: jk2635 !!----------------------------------------------------------------------2636 ! 2637 2638 zn1 = 1. /(jpk-1.)2639 zn2 = 1. - zn12640 2797 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2798 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2799 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2800 ! 2801 INTEGER :: jk ! dummy loop index 2802 REAL(wp) :: za1,za2,za3 ! local scalar 2803 REAL(wp) :: zn1,zn2 ! - - 2804 REAL(wp) :: za,zb,zx ! - - 2805 !!---------------------------------------------------------------------- 2806 ! 2807 zn1 = 1._wp / REAL( jpkm1, wp ) 2808 zn2 = 1._wp - zn1 2809 ! 2641 2810 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2642 2811 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2643 2812 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2644 2813 ! 2645 2814 za = pzb - za3*(pzs-za1)-za2 2646 2815 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2647 2816 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2648 2817 zx = 1.0_wp-za/2.0_wp-zb 2649 2818 ! 2650 2819 DO jk = 1, jpk 2651 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2652 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2653 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2820 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2821 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2822 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2654 2823 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2655 ENDDO 2656 2824 END DO 2657 2825 ! 2658 2826 END FUNCTION fgamma
Note: See TracChangeset
for help on using the changeset viewer.
