Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5836 r6140 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 … … 38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration 40 ! 40 41 USE in_out_manager ! I/O manager 41 42 USE iom ! I/O library … … 73 74 74 75 !! * Substitutions 75 # include "domzgr_substitute.h90"76 76 # include "vectopt_loop_substitute.h90" 77 77 !!---------------------------------------------------------------------- … … 102 102 INTEGER :: ios 103 103 ! 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 105 105 !!---------------------------------------------------------------------- 106 106 ! … … 120 120 WRITE(numout,*) 'dom_zgr : vertical coordinate' 121 121 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 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 127 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 127 128 ENDIF 129 130 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 131 129 132 ioptio = 0 ! Check Vertical coordinate options … … 155 158 ! 156 159 IF( nprint == 1 .AND. lwp ) THEN 157 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )160 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 158 161 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(:,:,:) )162 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 163 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 164 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 165 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 166 & ' w ', MINVAL( e3w_0(:,:,:) ) 164 167 165 168 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(:,:,:) )169 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 170 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 171 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 172 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 173 & ' w ', MAXVAL( e3w_0(:,:,:) ) 171 174 ENDIF 172 175 ! … … 379 382 !! - bathy : meter bathymetry (in meters) 380 383 !!---------------------------------------------------------------------- 381 INTEGER :: ji, jj, j l, jk ! dummy loop indices384 INTEGER :: ji, jj, jk ! dummy loop indices 382 385 INTEGER :: inum ! temporary logical unit 383 386 INTEGER :: ierror ! error flag … … 541 544 CALL iom_close( inum ) 542 545 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 546 547 ! set grounded point to 0 548 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 549 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 550 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 551 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 552 END WHERE 543 553 END IF 544 554 ! … … 578 588 ! 579 589 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 580 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin581 IF ( ln_isfcav ) THEN582 WHERE (bathy == risfdep)583 bathy = 0.0_wp ; risfdep = 0.0_wp584 END WHERE585 END IF586 ! end patch587 590 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 588 591 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 674 677 !! - update bathy : meter bathymetry (in meters) 675 678 !!---------------------------------------------------------------------- 676 !!677 679 INTEGER :: ji, jj, jl ! dummy loop indices 678 680 INTEGER :: icompt, ibtest, ikmax ! temporary integers 679 681 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 680 681 682 !!---------------------------------------------------------------------- 682 683 ! … … 775 776 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 776 777 ENDIF 777 778 IF( lwp .AND. nprint == 1 ) THEN ! control print779 WRITE(numout,*)780 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '781 WRITE(numout,*) ' ------------------'782 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )783 WRITE(numout,*)784 ENDIF785 778 ! 786 779 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 803 796 !! (min value = 1 over land) 804 797 !!---------------------------------------------------------------------- 805 !!806 798 INTEGER :: ji, jj ! dummy loop indices 807 799 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 835 827 END SUBROUTINE zgr_bot_level 836 828 837 SUBROUTINE zgr_top_level 838 !!---------------------------------------------------------------------- 839 !! *** ROUTINE zgr_bot_level *** 829 830 SUBROUTINE zgr_top_level 831 !!---------------------------------------------------------------------- 832 !! *** ROUTINE zgr_top_level *** 840 833 !! 841 834 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 847 840 !! (min value = 1) 848 841 !!---------------------------------------------------------------------- 849 !!850 842 INTEGER :: ji, jj ! dummy loop indices 851 843 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 881 873 END SUBROUTINE zgr_top_level 882 874 875 883 876 SUBROUTINE zgr_zco 884 877 !!---------------------------------------------------------------------- 885 878 !! *** ROUTINE zgr_zco *** 886 879 !! 887 !! ** Purpose : define the z-coordinate system880 !! ** Purpose : define the reference z-coordinate system 888 881 !! 889 882 !! ** Method : set 3D coord. arrays to reference 1D array … … 895 888 ! 896 889 DO jk = 1, jpk 897 gdept_0 898 gdepw_0 899 gde p3w_0(:,:,jk) = gdepw_1d(jk)900 e3t_0 901 e3u_0 902 e3v_0 903 e3f_0 904 e3w_0 905 e3uw_0 906 e3vw_0 890 gdept_0(:,:,jk) = gdept_1d(jk) 891 gdepw_0(:,:,jk) = gdepw_1d(jk) 892 gde3w_0(:,:,jk) = gdepw_1d(jk) 893 e3t_0 (:,:,jk) = e3t_1d (jk) 894 e3u_0 (:,:,jk) = e3t_1d (jk) 895 e3v_0 (:,:,jk) = e3t_1d (jk) 896 e3f_0 (:,:,jk) = e3t_1d (jk) 897 e3w_0 (:,:,jk) = e3w_1d (jk) 898 e3uw_0 (:,:,jk) = e3w_1d (jk) 899 e3vw_0 (:,:,jk) = e3w_1d (jk) 907 900 END DO 908 901 ! … … 917 910 !! 918 911 !! ** Purpose : the depth and vertical scale factor in partial step 919 !! z-coordinate case912 !! reference z-coordinate case 920 913 !! 921 914 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 957 950 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 958 951 !!---------------------------------------------------------------------- 959 !!960 952 INTEGER :: ji, jj, jk ! dummy loop indices 961 953 INTEGER :: ik, it, ikb, ikt ! temporary integers 962 LOGICAL :: ll_print ! Allow control print for debugging963 954 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 964 955 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 965 REAL(wp) :: zmax ! Maximum depth966 956 REAL(wp) :: zdiff ! temporary scalar 967 REAL(wp) :: z refdep! temporary scalar957 REAL(wp) :: zmax ! temporary scalar 968 958 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 969 959 !!--------------------------------------------------------------------- … … 971 961 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 972 962 ! 973 CALL wrk_alloc( jpi, jpj, jpk,zprt )963 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 974 964 ! 975 965 IF(lwp) WRITE(numout,*) … … 977 967 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 978 968 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 979 980 ll_print = .FALSE. ! Local variable for debugging981 982 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth983 WRITE(numout,*)984 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'985 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )986 ENDIF987 988 969 989 970 ! bathymetry in level (from bathy_meter) … … 1004 985 END DO 1005 986 1006 IF ( ln_isfcav ) CALL zgr_isf1007 1008 987 ! Scale factors and depth at T- and W-points 1009 988 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1013 992 e3w_0 (:,:,jk) = e3w_1d (jk) 1014 993 END DO 994 995 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 996 IF ( ln_isfcav ) CALL zgr_isf 997 998 ! Scale factors and depth at T- and W-points 999 IF ( .NOT. ln_isfcav ) THEN 1000 DO jj = 1, jpj 1001 DO ji = 1, jpi 1002 ik = mbathy(ji,jj) 1003 IF( ik > 0 ) THEN ! ocean point only 1004 ! max ocean level case 1005 IF( ik == jpkm1 ) THEN 1006 zdepwp = bathy(ji,jj) 1007 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1008 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1009 e3t_0(ji,jj,ik ) = ze3tp 1010 e3t_0(ji,jj,ik+1) = ze3tp 1011 e3w_0(ji,jj,ik ) = ze3wp 1012 e3w_0(ji,jj,ik+1) = ze3tp 1013 gdepw_0(ji,jj,ik+1) = zdepwp 1014 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1015 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1016 ! 1017 ELSE ! standard case 1018 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1019 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1020 ENDIF 1021 !gm Bug? check the gdepw_1d 1022 ! ... on ik 1023 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1024 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1025 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1026 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1027 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1028 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1029 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1030 ! ... on ik+1 1031 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1032 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1033 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1034 ENDIF 1035 ENDIF 1036 END DO 1037 END DO 1038 ! 1039 it = 0 1040 DO jj = 1, jpj 1041 DO ji = 1, jpi 1042 ik = mbathy(ji,jj) 1043 IF( ik > 0 ) THEN ! ocean point only 1044 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1045 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1046 ! test 1047 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1048 IF( zdiff <= 0._wp .AND. lwp ) THEN 1049 it = it + 1 1050 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1051 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1052 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1053 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1054 ENDIF 1055 ENDIF 1056 END DO 1057 END DO 1058 END IF 1059 ! 1060 ! Scale factors and depth at U-, V-, UW and VW-points 1061 DO jk = 1, jpk ! initialisation to z-scale factors 1062 e3u_0 (:,:,jk) = e3t_1d(jk) 1063 e3v_0 (:,:,jk) = e3t_1d(jk) 1064 e3uw_0(:,:,jk) = e3w_1d(jk) 1065 e3vw_0(:,:,jk) = e3w_1d(jk) 1066 END DO 1067 1068 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1069 DO jj = 1, jpjm1 1070 DO ji = 1, fs_jpim1 ! vector opt. 1071 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1072 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1073 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1074 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1075 END DO 1076 END DO 1077 END DO 1078 IF ( ln_isfcav ) THEN 1079 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1080 DO jj = 2, jpjm1 1081 DO ji = 2, fs_jpim1 ! vector opt. 1082 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1083 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1084 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1085 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1086 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1087 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1088 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1089 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1090 END DO 1091 END DO 1092 END IF 1093 1094 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1095 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1096 ! 1097 1098 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1099 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1100 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1101 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1102 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1103 END DO 1104 1105 ! Scale factor at F-point 1106 DO jk = 1, jpk ! initialisation to z-scale factors 1107 e3f_0(:,:,jk) = e3t_1d(jk) 1108 END DO 1109 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1110 DO jj = 1, jpjm1 1111 DO ji = 1, fs_jpim1 ! vector opt. 1112 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1113 END DO 1114 END DO 1115 END DO 1116 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1117 ! 1118 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1119 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1120 END DO 1121 !!gm bug ? : must be a do loop with mj0,mj1 1015 1122 ! 1123 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1124 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1125 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1126 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1127 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1128 1129 ! Control of the sign 1130 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1131 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1132 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1133 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1134 1135 ! Compute gde3w_0 (vertical sum of e3w) 1136 IF ( ln_isfcav ) THEN ! if cavity 1137 WHERE( misfdep == 0 ) misfdep = 1 1138 DO jj = 1,jpj 1139 DO ji = 1,jpi 1140 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1141 DO jk = 2, misfdep(ji,jj) 1142 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1143 END DO 1144 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)) 1145 DO jk = misfdep(ji,jj) + 1, jpk 1146 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1147 END DO 1148 END DO 1149 END DO 1150 ELSE ! no cavity 1151 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1152 DO jk = 2, jpk 1153 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1154 END DO 1155 END IF 1156 ! 1157 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1158 ! 1159 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1160 ! 1161 END SUBROUTINE zgr_zps 1162 1163 1164 SUBROUTINE zgr_isf 1165 !!---------------------------------------------------------------------- 1166 !! *** ROUTINE zgr_isf *** 1167 !! 1168 !! ** Purpose : check the bathymetry in levels 1169 !! 1170 !! ** Method : THe water column have to contained at least 2 cells 1171 !! Bathymetry and isfdraft are modified (dig/close) to respect 1172 !! this criterion. 1173 !! 1174 !! ** Action : - test compatibility between isfdraft and bathy 1175 !! - bathy and isfdraft are modified 1176 !!---------------------------------------------------------------------- 1177 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1178 INTEGER :: ik, it ! temporary integers 1179 INTEGER :: icompt, ibtest ! (ISF) 1180 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1181 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1182 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1183 REAL(wp) :: zmax ! Maximum and minimum depth 1184 REAL(wp) :: zbathydiff ! isf temporary scalar 1185 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1186 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1187 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1188 REAL(wp) :: zdiff ! temporary scalar 1189 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1190 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1191 !!--------------------------------------------------------------------- 1192 ! 1193 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1194 ! 1195 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1196 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1197 1198 ! (ISF) compute misfdep 1199 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1200 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1201 END WHERE 1202 1203 ! Compute misfdep for ocean points (i.e. first wet level) 1204 ! find the first ocean level such that the first level thickness 1205 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1206 ! e3t_0 is the reference level thickness 1207 DO jk = 2, jpkm1 1208 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1209 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1210 END DO 1211 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1212 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1213 END WHERE 1214 1215 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1216 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1217 misfdep = 0; risfdep = 0.0_wp; 1218 mbathy = 0; bathy = 0.0_wp; 1219 END WHERE 1220 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1221 misfdep = 0; risfdep = 0.0_wp; 1222 mbathy = 0; bathy = 0.0_wp; 1223 END WHERE 1224 1225 ! 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 1226 icompt = 0 1227 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1228 DO jl = 1, 10 1229 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1230 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1231 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1232 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1233 END WHERE 1234 WHERE (mbathy(:,:) <= 0) 1235 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1236 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1237 END WHERE 1238 IF( lk_mpp ) THEN 1239 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1240 CALL lbc_lnk( zbathy, 'T', 1. ) 1241 misfdep(:,:) = INT( zbathy(:,:) ) 1242 1243 CALL lbc_lnk( risfdep,'T', 1. ) 1244 CALL lbc_lnk( bathy, 'T', 1. ) 1245 1246 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1247 CALL lbc_lnk( zbathy, 'T', 1. ) 1248 mbathy(:,:) = INT( zbathy(:,:) ) 1249 ENDIF 1250 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1251 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1252 misfdep(jpi,:) = misfdep( 2 ,:) 1253 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1254 mbathy(jpi,:) = mbathy( 2 ,:) 1255 ENDIF 1256 1257 ! split last cell if possible (only where water column is 2 cell or less) 1258 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1259 IF ( .NOT. ln_iscpl) THEN 1260 DO jk = jpkm1, 1, -1 1261 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1262 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1263 mbathy(:,:) = jk 1264 bathy(:,:) = zmax 1265 END WHERE 1266 END DO 1267 END IF 1268 1269 ! split top cell if possible (only where water column is 2 cell or less) 1270 DO jk = 2, jpkm1 1271 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1272 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1273 misfdep(:,:) = jk 1274 risfdep(:,:) = zmax 1275 END WHERE 1276 END DO 1277 1278 1279 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1280 DO jj = 1, jpj 1281 DO ji = 1, jpi 1282 ! find the minimum change option: 1283 ! test bathy 1284 IF (risfdep(ji,jj) > 1) THEN 1285 IF ( .NOT. ln_iscpl ) THEN 1286 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1287 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1288 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1289 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1290 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1291 IF (zbathydiff <= zrisfdepdiff) THEN 1292 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1293 mbathy(ji,jj)= mbathy(ji,jj) + 1 1294 ELSE 1295 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1296 misfdep(ji,jj) = misfdep(ji,jj) - 1 1297 END IF 1298 ENDIF 1299 ELSE 1300 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1301 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1302 misfdep(ji,jj) = misfdep(ji,jj) - 1 1303 END IF 1304 END IF 1305 END IF 1306 END DO 1307 END DO 1308 1309 ! At least 2 levels for water thickness at T, U, and V point. 1310 DO jj = 1, jpj 1311 DO ji = 1, jpi 1312 ! find the minimum change option: 1313 ! test bathy 1314 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1315 IF ( .NOT. ln_iscpl ) THEN 1316 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1317 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1318 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1319 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1320 IF (zbathydiff <= zrisfdepdiff) THEN 1321 mbathy(ji,jj) = mbathy(ji,jj) + 1 1322 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1323 ELSE 1324 misfdep(ji,jj)= misfdep(ji,jj) - 1 1325 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1326 END IF 1327 ELSE 1328 misfdep(ji,jj)= misfdep(ji,jj) - 1 1329 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1330 END IF 1331 ENDIF 1332 END DO 1333 END DO 1334 1335 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1336 DO jj = 1, jpjm1 1337 DO ji = 1, jpim1 1338 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1339 IF ( .NOT. ln_iscpl ) THEN 1340 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1341 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1342 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1343 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1344 IF (zbathydiff <= zrisfdepdiff) THEN 1345 mbathy(ji,jj) = mbathy(ji,jj) + 1 1346 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1347 ELSE 1348 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1349 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1350 END IF 1351 ELSE 1352 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1353 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1354 END IF 1355 ENDIF 1356 END DO 1357 END DO 1358 1359 IF( lk_mpp ) THEN 1360 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1361 CALL lbc_lnk( zbathy, 'T', 1. ) 1362 misfdep(:,:) = INT( zbathy(:,:) ) 1363 1364 CALL lbc_lnk( risfdep,'T', 1. ) 1365 CALL lbc_lnk( bathy, 'T', 1. ) 1366 1367 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1368 CALL lbc_lnk( zbathy, 'T', 1. ) 1369 mbathy(:,:) = INT( zbathy(:,:) ) 1370 ENDIF 1371 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1372 DO jj = 1, jpjm1 1373 DO ji = 1, jpim1 1374 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1375 IF ( .NOT. ln_iscpl ) THEN 1376 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1377 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1378 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1379 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1380 IF (zbathydiff <= zrisfdepdiff) THEN 1381 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1382 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1383 ELSE 1384 misfdep(ji,jj) = misfdep(ji,jj) - 1 1385 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1386 END IF 1387 ELSE 1388 misfdep(ji,jj) = misfdep(ji,jj) - 1 1389 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1390 END IF 1391 ENDIF 1392 END DO 1393 END DO 1394 1395 1396 IF( lk_mpp ) THEN 1397 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1398 CALL lbc_lnk( zbathy, 'T', 1. ) 1399 misfdep(:,:) = INT( zbathy(:,:) ) 1400 1401 CALL lbc_lnk( risfdep,'T', 1. ) 1402 CALL lbc_lnk( bathy, 'T', 1. ) 1403 1404 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1405 CALL lbc_lnk( zbathy, 'T', 1. ) 1406 mbathy(:,:) = INT( zbathy(:,:) ) 1407 ENDIF 1408 1409 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1410 DO jj = 1, jpjm1 1411 DO ji = 1, jpim1 1412 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1413 IF ( .NOT. ln_iscpl ) THEN 1414 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+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1417 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1418 IF (zbathydiff <= zrisfdepdiff) THEN 1419 mbathy(ji,jj) = mbathy(ji,jj) + 1 1420 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1421 ELSE 1422 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1423 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1424 END IF 1425 ELSE 1426 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1427 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1428 ENDIF 1429 ENDIF 1430 ENDDO 1431 ENDDO 1432 1433 IF( lk_mpp ) THEN 1434 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1435 CALL lbc_lnk( zbathy, 'T', 1. ) 1436 misfdep(:,:) = INT( zbathy(:,:) ) 1437 1438 CALL lbc_lnk( risfdep,'T', 1. ) 1439 CALL lbc_lnk( bathy, 'T', 1. ) 1440 1441 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1442 CALL lbc_lnk( zbathy, 'T', 1. ) 1443 mbathy(:,:) = INT( zbathy(:,:) ) 1444 ENDIF 1445 1446 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1447 DO jj = 1, jpjm1 1448 DO ji = 1, jpim1 1449 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1450 IF ( .NOT. ln_iscpl ) THEN 1451 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1452 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1453 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1454 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1455 IF (zbathydiff <= zrisfdepdiff) THEN 1456 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1457 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1458 ELSE 1459 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1460 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1461 END IF 1462 ELSE 1463 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1464 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1465 ENDIF 1466 ENDIF 1467 ENDDO 1468 ENDDO 1469 1470 IF( lk_mpp ) THEN 1471 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1472 CALL lbc_lnk( zbathy, 'T', 1. ) 1473 misfdep(:,:) = INT( zbathy(:,:) ) 1474 1475 CALL lbc_lnk( risfdep,'T', 1. ) 1476 CALL lbc_lnk( bathy, 'T', 1. ) 1477 1478 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1479 CALL lbc_lnk( zbathy, 'T', 1. ) 1480 mbathy(:,:) = INT( zbathy(:,:) ) 1481 ENDIF 1482 END DO 1483 ! end dig bathy/ice shelf to be compatible 1484 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1485 DO jl = 1,20 1486 1487 ! remove single point "bay" on isf coast line in the ice shelf draft' 1488 DO jk = 2, jpk 1489 WHERE (misfdep==0) misfdep=jpk 1490 zmask=0._wp 1491 WHERE (misfdep <= jk) zmask=1 1492 DO jj = 2, jpjm1 1493 DO ji = 2, jpim1 1494 IF (misfdep(ji,jj) == jk) THEN 1495 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1496 IF (ibtest <= 1) THEN 1497 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1498 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1499 END IF 1500 END IF 1501 END DO 1502 END DO 1503 END DO 1504 WHERE (misfdep==jpk) 1505 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1506 END WHERE 1507 IF( lk_mpp ) THEN 1508 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1509 CALL lbc_lnk( zbathy, 'T', 1. ) 1510 misfdep(:,:) = INT( zbathy(:,:) ) 1511 1512 CALL lbc_lnk( risfdep,'T', 1. ) 1513 CALL lbc_lnk( bathy, 'T', 1. ) 1514 1515 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1516 CALL lbc_lnk( zbathy, 'T', 1. ) 1517 mbathy(:,:) = INT( zbathy(:,:) ) 1518 ENDIF 1519 1520 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1521 DO jk = jpk,1,-1 1522 zmask=0._wp 1523 WHERE (mbathy >= jk ) zmask=1 1524 DO jj = 2, jpjm1 1525 DO ji = 2, jpim1 1526 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1527 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1528 IF (ibtest <= 1) THEN 1529 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1530 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1531 END IF 1532 END IF 1533 END DO 1534 END DO 1535 END DO 1536 WHERE (mbathy==0) 1537 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1538 END WHERE 1539 IF( lk_mpp ) THEN 1540 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1541 CALL lbc_lnk( zbathy, 'T', 1. ) 1542 misfdep(:,:) = INT( zbathy(:,:) ) 1543 1544 CALL lbc_lnk( risfdep,'T', 1. ) 1545 CALL lbc_lnk( bathy, 'T', 1. ) 1546 1547 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1548 CALL lbc_lnk( zbathy, 'T', 1. ) 1549 mbathy(:,:) = INT( zbathy(:,:) ) 1550 ENDIF 1551 1552 ! fill hole in ice shelf 1553 zmisfdep = misfdep 1554 zrisfdep = risfdep 1555 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1556 DO jj = 2, jpjm1 1557 DO ji = 2, jpim1 1558 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1559 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1560 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1561 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1562 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1563 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1564 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1565 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1566 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1567 END IF 1568 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1569 misfdep(ji,jj) = ibtest 1570 risfdep(ji,jj) = gdepw_1d(ibtest) 1571 ENDIF 1572 ENDDO 1573 ENDDO 1574 1575 IF( lk_mpp ) THEN 1576 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1577 CALL lbc_lnk( zbathy, 'T', 1. ) 1578 misfdep(:,:) = INT( zbathy(:,:) ) 1579 1580 CALL lbc_lnk( risfdep, 'T', 1. ) 1581 CALL lbc_lnk( bathy, 'T', 1. ) 1582 1583 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1584 CALL lbc_lnk( zbathy, 'T', 1. ) 1585 mbathy(:,:) = INT( zbathy(:,:) ) 1586 ENDIF 1587 ! 1588 !! fill hole in bathymetry 1589 zmbathy (:,:)=mbathy (:,:) 1590 DO jj = 2, jpjm1 1591 DO ji = 2, jpim1 1592 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1593 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1594 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1595 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1596 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1597 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1598 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1599 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1600 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1601 END IF 1602 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1603 mbathy(ji,jj) = ibtest 1604 bathy(ji,jj) = gdepw_1d(ibtest+1) 1605 ENDIF 1606 END DO 1607 END DO 1608 IF( lk_mpp ) THEN 1609 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1610 CALL lbc_lnk( zbathy, 'T', 1. ) 1611 misfdep(:,:) = INT( zbathy(:,:) ) 1612 1613 CALL lbc_lnk( risfdep, 'T', 1. ) 1614 CALL lbc_lnk( bathy, 'T', 1. ) 1615 1616 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1617 CALL lbc_lnk( zbathy, 'T', 1. ) 1618 mbathy(:,:) = INT( zbathy(:,:) ) 1619 ENDIF 1620 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1621 DO jj = 1, jpjm1 1622 DO ji = 1, jpim1 1623 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1624 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1625 END IF 1626 END DO 1627 END DO 1628 IF( lk_mpp ) THEN 1629 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1630 CALL lbc_lnk( zbathy, 'T', 1. ) 1631 misfdep(:,:) = INT( zbathy(:,:) ) 1632 1633 CALL lbc_lnk( risfdep, 'T', 1. ) 1634 CALL lbc_lnk( bathy, 'T', 1. ) 1635 1636 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1637 CALL lbc_lnk( zbathy, 'T', 1. ) 1638 mbathy(:,:) = INT( zbathy(:,:) ) 1639 ENDIF 1640 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1641 DO jj = 1, jpjm1 1642 DO ji = 1, jpim1 1643 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1644 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1645 END IF 1646 END DO 1647 END DO 1648 IF( lk_mpp ) THEN 1649 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1650 CALL lbc_lnk( zbathy, 'T', 1. ) 1651 misfdep(:,:) = INT( zbathy(:,:) ) 1652 1653 CALL lbc_lnk( risfdep,'T', 1. ) 1654 CALL lbc_lnk( bathy, 'T', 1. ) 1655 1656 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1657 CALL lbc_lnk( zbathy, 'T', 1. ) 1658 mbathy(:,:) = INT( zbathy(:,:) ) 1659 ENDIF 1660 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1661 DO jj = 1, jpjm1 1662 DO ji = 1, jpi 1663 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1664 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1665 END IF 1666 END DO 1667 END DO 1668 IF( lk_mpp ) THEN 1669 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1670 CALL lbc_lnk( zbathy, 'T', 1. ) 1671 misfdep(:,:) = INT( zbathy(:,:) ) 1672 1673 CALL lbc_lnk( risfdep,'T', 1. ) 1674 CALL lbc_lnk( bathy, 'T', 1. ) 1675 1676 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1677 CALL lbc_lnk( zbathy, 'T', 1. ) 1678 mbathy(:,:) = INT( zbathy(:,:) ) 1679 ENDIF 1680 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1681 DO jj = 1, jpjm1 1682 DO ji = 1, jpi 1683 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1684 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1685 END IF 1686 END DO 1687 END DO 1688 IF( lk_mpp ) THEN 1689 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1690 CALL lbc_lnk( zbathy, 'T', 1. ) 1691 misfdep(:,:) = INT( zbathy(:,:) ) 1692 1693 CALL lbc_lnk( risfdep,'T', 1. ) 1694 CALL lbc_lnk( bathy, 'T', 1. ) 1695 1696 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1697 CALL lbc_lnk( zbathy, 'T', 1. ) 1698 mbathy(:,:) = INT( zbathy(:,:) ) 1699 ENDIF 1700 ! if not compatible after all check, mask T 1701 DO jj = 1, jpj 1702 DO ji = 1, jpi 1703 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1704 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1705 END IF 1706 END DO 1707 END DO 1708 1709 WHERE (mbathy(:,:) == 1) 1710 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1711 END WHERE 1712 END DO 1713 ! end check compatibility ice shelf/bathy 1714 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1715 WHERE (risfdep(:,:) <= 10._wp) 1716 misfdep = 1; risfdep = 0.0_wp; 1717 END WHERE 1718 1719 IF( icompt == 0 ) THEN 1720 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1721 ELSE 1722 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1723 ENDIF 1724 1725 ! compute scale factor and depth at T- and W- points 1016 1726 DO jj = 1, jpj 1017 1727 DO ji = 1, jpi … … 1035 1745 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1036 1746 ENDIF 1747 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1037 1748 !gm Bug? check the gdepw_1d 1038 1749 ! ... on ik … … 1040 1751 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1041 1752 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1042 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1043 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1044 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1045 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1753 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1754 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1046 1755 ! ... on ik+1 1047 1756 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1048 1757 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1049 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1050 1758 ENDIF 1051 1759 ENDIF … … 1073 1781 END DO 1074 1782 ! 1075 IF ( ln_isfcav ) THEN1076 1783 ! (ISF) Definition of e3t, u, v, w for ISF case 1077 1078 1079 1080 1081 1082 1784 DO jj = 1, jpj 1785 DO ji = 1, jpi 1786 ik = misfdep(ji,jj) 1787 IF( ik > 1 ) THEN ! ice shelf point only 1788 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1789 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1083 1790 !gm Bug? check the gdepw_0 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 e3w_0 (ji,jj,ik ) =2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1096 1791 ! ... on ik 1792 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1793 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1794 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1795 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1796 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1797 1798 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1799 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1800 ENDIF 1801 ! ... on ik / ik-1 1802 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1803 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1097 1804 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1098 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1805 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1806 ENDIF 1807 END DO 1808 END DO 1809 1810 it = 0 1811 DO jj = 1, jpj 1812 DO ji = 1, jpi 1813 ik = misfdep(ji,jj) 1814 IF( ik > 1 ) THEN ! ice shelf point only 1815 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1816 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1817 ! test 1818 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1819 IF( zdiff <= 0. .AND. lwp ) THEN 1820 it = it + 1 1821 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1822 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1823 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1824 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1099 1825 ENDIF 1100 END DO1826 ENDIF 1101 1827 END DO 1102 !1103 it = 01104 DO jj = 1, jpj1105 DO ji = 1, jpi1106 ik = misfdep(ji,jj)1107 IF( ik > 1 ) THEN ! ice shelf point only1108 e3tp (ji,jj) = e3t_0(ji,jj,ik )1109 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1110 ! test1111 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1112 IF( zdiff <= 0. .AND. lwp ) THEN1113 it = it + 11114 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1115 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1116 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1117 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1118 ENDIF1119 ENDIF1120 END DO1121 END DO1122 END IF1123 ! END (ISF)1124 1125 ! Scale factors and depth at U-, V-, UW and VW-points1126 DO jk = 1, jpk ! initialisation to z-scale factors1127 e3u_0 (:,:,jk) = e3t_1d(jk)1128 e3v_0 (:,:,jk) = e3t_1d(jk)1129 e3uw_0(:,:,jk) = e3w_1d(jk)1130 e3vw_0(:,:,jk) = e3w_1d(jk)1131 END DO1132 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1133 DO jj = 1, jpjm11134 DO ji = 1, fs_jpim1 ! vector opt.1135 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1136 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1137 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1138 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1139 END DO1140 END DO1141 END DO1142 IF ( ln_isfcav ) THEN1143 ! (ISF) define e3uw (adapted for 2 cells in the water column)1144 DO jj = 2, jpjm11145 DO ji = 2, fs_jpim1 ! vector opt.1146 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1147 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1148 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1149 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1150 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1151 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1152 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1153 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1154 END DO1155 END DO1156 END IF1157 1158 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1159 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1160 !1161 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1162 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1163 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1164 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1165 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1166 END DO1167 1168 ! Scale factor at F-point1169 DO jk = 1, jpk ! initialisation to z-scale factors1170 e3f_0(:,:,jk) = e3t_1d(jk)1171 END DO1172 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1173 DO jj = 1, jpjm11174 DO ji = 1, fs_jpim1 ! vector opt.1175 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1176 END DO1177 END DO1178 END DO1179 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1180 !1181 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1182 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1183 END DO1184 !!gm bug ? : must be a do loop with mj0,mj11185 !1186 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21187 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1188 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1189 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1190 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1191 1192 ! Control of the sign1193 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1194 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1195 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1196 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1197 1198 ! Compute gdep3w_0 (vertical sum of e3w)1199 IF ( ln_isfcav ) THEN ! if cavity1200 WHERE (misfdep == 0) misfdep = 11201 DO jj = 1,jpj1202 DO ji = 1,jpi1203 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1204 DO jk = 2, misfdep(ji,jj)1205 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1206 END DO1207 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))1208 DO jk = misfdep(ji,jj) + 1, jpk1209 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1210 END DO1211 END DO1212 END DO1213 ELSE ! no cavity1214 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1215 DO jk = 2, jpk1216 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1217 END DO1218 END IF1219 ! ! ================= !1220 IF(lwp .AND. ll_print) THEN ! Control print !1221 ! ! ================= !1222 DO jj = 1,jpj1223 DO ji = 1, jpi1224 ik = MAX( mbathy(ji,jj), 1 )1225 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1226 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1227 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1228 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1229 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1230 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1231 END DO1232 END DO1233 WRITE(numout,*)1234 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1235 WRITE(numout,*)1236 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1237 WRITE(numout,*)1238 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1239 WRITE(numout,*)1240 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1241 WRITE(numout,*)1242 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1243 WRITE(numout,*)1244 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1245 ENDIF1246 !1247 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1248 !1249 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1250 !1251 END SUBROUTINE zgr_zps1252 1253 SUBROUTINE zgr_isf1254 !!----------------------------------------------------------------------1255 !! *** ROUTINE zgr_isf ***1256 !!1257 !! ** Purpose : check the bathymetry in levels1258 !!1259 !! ** Method : THe water column have to contained at least 2 cells1260 !! Bathymetry and isfdraft are modified (dig/close) to respect1261 !! this criterion.1262 !!1263 !!1264 !! ** Action : - test compatibility between isfdraft and bathy1265 !! - bathy and isfdraft are modified1266 !!----------------------------------------------------------------------1267 !!1268 INTEGER :: ji, jj, jk, jl ! dummy loop indices1269 INTEGER :: ik, it ! temporary integers1270 INTEGER :: id, jd, nprocd1271 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1272 LOGICAL :: ll_print ! Allow control print for debugging1273 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1274 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1275 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1276 REAL(wp) :: zdiff ! temporary scalar1277 REAL(wp) :: zrefdep ! temporary scalar1278 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1279 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1280 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1281 !!---------------------------------------------------------------------1282 !1283 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1284 !1285 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)1286 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )1287 1288 1289 ! (ISF) compute misfdep1290 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11291 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1292 END WHERE1293 1294 ! Compute misfdep for ocean points (i.e. first wet level)1295 ! find the first ocean level such that the first level thickness1296 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1297 ! e3t_0 is the reference level thickness1298 DO jk = 2, jpkm11299 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1300 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11301 1828 END DO 1302 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1303 risfdep(:,:) = 0. ; misfdep(:,:) = 11304 END WHERE1305 1306 ! 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 situation1307 icompt = 01308 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1309 DO jl = 1, 101310 WHERE (bathy(:,:) .EQ. risfdep(:,:) )1311 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1312 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1313 END WHERE1314 WHERE (mbathy(:,:) <= 0)1315 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1316 mbathy (:,:) = 0; bathy (:,:) = 0._wp1317 ENDWHERE1318 IF( lk_mpp ) THEN1319 zbathy(:,:) = FLOAT( misfdep(:,:) )1320 CALL lbc_lnk( zbathy, 'T', 1. )1321 misfdep(:,:) = INT( zbathy(:,:) )1322 CALL lbc_lnk( risfdep, 'T', 1. )1323 CALL lbc_lnk( bathy, 'T', 1. )1324 zbathy(:,:) = FLOAT( mbathy(:,:) )1325 CALL lbc_lnk( zbathy, 'T', 1. )1326 mbathy(:,:) = INT( zbathy(:,:) )1327 ENDIF1328 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1329 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1330 misfdep(jpi,:) = misfdep( 2 ,:)1331 ENDIF1332 1333 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1334 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1335 mbathy(jpi,:) = mbathy( 2 ,:)1336 ENDIF1337 1338 ! split last cell if possible (only where water column is 2 cell or less)1339 DO jk = jpkm1, 1, -11340 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1341 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1342 mbathy(:,:) = jk1343 bathy(:,:) = zmax1344 END WHERE1345 END DO1346 1347 ! split top cell if possible (only where water column is 2 cell or less)1348 DO jk = 2, jpkm11349 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1350 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1351 misfdep(:,:) = jk1352 risfdep(:,:) = zmax1353 END WHERE1354 END DO1355 1356 1357 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1358 DO jj = 1, jpj1359 DO ji = 1, jpi1360 ! find the minimum change option:1361 ! test bathy1362 IF (risfdep(ji,jj) .GT. 1) THEN1363 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1364 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1365 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1366 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1367 1368 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN1369 IF (zbathydiff .LE. zrisfdepdiff) THEN1370 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1371 mbathy(ji,jj)= mbathy(ji,jj) + 11372 ELSE1373 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1374 misfdep(ji,jj) = misfdep(ji,jj) - 11375 END IF1376 END IF1377 END IF1378 END DO1379 END DO1380 1381 ! At least 2 levels for water thickness at T, U, and V point.1382 DO jj = 1, jpj1383 DO ji = 1, jpi1384 ! find the minimum change option:1385 ! test bathy1386 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1387 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1388 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1389 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1390 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1391 IF (zbathydiff .LE. zrisfdepdiff) THEN1392 mbathy(ji,jj) = mbathy(ji,jj) + 11393 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1394 ELSE1395 misfdep(ji,jj)= misfdep(ji,jj) - 11396 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1397 END IF1398 ENDIF1399 END DO1400 END DO1401 1402 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1403 DO jj = 1, jpjm11404 DO ji = 1, jpim11405 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1406 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1407 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1408 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1409 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1410 IF (zbathydiff .LE. zrisfdepdiff) THEN1411 mbathy(ji,jj) = mbathy(ji,jj) + 11412 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1413 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1414 ELSE1415 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11416 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1417 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1418 END IF1419 ENDIF1420 END DO1421 END DO1422 1423 IF( lk_mpp ) THEN1424 zbathy(:,:) = FLOAT( misfdep(:,:) )1425 CALL lbc_lnk( zbathy, 'T', 1. )1426 misfdep(:,:) = INT( zbathy(:,:) )1427 CALL lbc_lnk( risfdep, 'T', 1. )1428 CALL lbc_lnk( bathy, 'T', 1. )1429 zbathy(:,:) = FLOAT( mbathy(:,:) )1430 CALL lbc_lnk( zbathy, 'T', 1. )1431 mbathy(:,:) = INT( zbathy(:,:) )1432 ENDIF1433 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1434 DO jj = 1, jpjm11435 DO ji = 1, jpim11436 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1437 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1438 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1439 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1440 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1441 IF (zbathydiff .LE. zrisfdepdiff) THEN1442 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11443 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1444 ELSE1445 misfdep(ji,jj) = misfdep(ji,jj) - 11446 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1447 END IF1448 ENDIF1449 END DO1450 END DO1451 1452 1453 IF( lk_mpp ) THEN1454 zbathy(:,:) = FLOAT( misfdep(:,:) )1455 CALL lbc_lnk( zbathy, 'T', 1. )1456 misfdep(:,:) = INT( zbathy(:,:) )1457 CALL lbc_lnk( risfdep, 'T', 1. )1458 CALL lbc_lnk( bathy, 'T', 1. )1459 zbathy(:,:) = FLOAT( mbathy(:,:) )1460 CALL lbc_lnk( zbathy, 'T', 1. )1461 mbathy(:,:) = INT( zbathy(:,:) )1462 ENDIF1463 1464 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1465 DO jj = 1, jpjm11466 DO ji = 1, jpim11467 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1468 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1469 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1470 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1471 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1472 IF (zbathydiff .LE. zrisfdepdiff) THEN1473 mbathy(ji,jj) = mbathy(ji,jj) + 11474 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1475 ELSE1476 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11477 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1478 END IF1479 ENDIF1480 ENDDO1481 ENDDO1482 1483 IF( lk_mpp ) THEN1484 zbathy(:,:) = FLOAT( misfdep(:,:) )1485 CALL lbc_lnk( zbathy, 'T', 1. )1486 misfdep(:,:) = INT( zbathy(:,:) )1487 CALL lbc_lnk( risfdep, 'T', 1. )1488 CALL lbc_lnk( bathy, 'T', 1. )1489 zbathy(:,:) = FLOAT( mbathy(:,:) )1490 CALL lbc_lnk( zbathy, 'T', 1. )1491 mbathy(:,:) = INT( zbathy(:,:) )1492 ENDIF1493 1494 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1495 DO jj = 1, jpjm11496 DO ji = 1, jpim11497 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1498 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1499 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1500 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1501 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1502 IF (zbathydiff .LE. zrisfdepdiff) THEN1503 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11504 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1505 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1506 ELSE1507 misfdep(ji,jj) = misfdep(ji ,jj) - 11508 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1509 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1510 END IF1511 ENDIF1512 ENDDO1513 ENDDO1514 1515 IF( lk_mpp ) THEN1516 zbathy(:,:) = FLOAT( misfdep(:,:) )1517 CALL lbc_lnk( zbathy, 'T', 1. )1518 misfdep(:,:) = INT( zbathy(:,:) )1519 CALL lbc_lnk( risfdep, 'T', 1. )1520 CALL lbc_lnk( bathy, 'T', 1. )1521 zbathy(:,:) = FLOAT( mbathy(:,:) )1522 CALL lbc_lnk( zbathy, 'T', 1. )1523 mbathy(:,:) = INT( zbathy(:,:) )1524 ENDIF1525 END DO1526 ! end dig bathy/ice shelf to be compatible1527 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1528 DO jl = 1,201529 1530 ! remove single point "bay" on isf coast line in the ice shelf draft'1531 DO jk = 2, jpk1532 WHERE (misfdep==0) misfdep=jpk1533 zmask=01534 WHERE (misfdep .LE. jk) zmask=11535 DO jj = 2, jpjm11536 DO ji = 2, jpim11537 IF (misfdep(ji,jj) .EQ. jk) THEN1538 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1539 IF (ibtest .LE. 1) THEN1540 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11541 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1542 END IF1543 END IF1544 END DO1545 END DO1546 END DO1547 WHERE (misfdep==jpk)1548 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1549 END WHERE1550 IF( lk_mpp ) THEN1551 zbathy(:,:) = FLOAT( misfdep(:,:) )1552 CALL lbc_lnk( zbathy, 'T', 1. )1553 misfdep(:,:) = INT( zbathy(:,:) )1554 CALL lbc_lnk( risfdep, 'T', 1. )1555 CALL lbc_lnk( bathy, 'T', 1. )1556 zbathy(:,:) = FLOAT( mbathy(:,:) )1557 CALL lbc_lnk( zbathy, 'T', 1. )1558 mbathy(:,:) = INT( zbathy(:,:) )1559 ENDIF1560 1561 ! remove single point "bay" on bathy coast line beneath an ice shelf'1562 DO jk = jpk,1,-11563 zmask=01564 WHERE (mbathy .GE. jk ) zmask=11565 DO jj = 2, jpjm11566 DO ji = 2, jpim11567 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1568 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1569 IF (ibtest .LE. 1) THEN1570 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11571 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01572 END IF1573 END IF1574 END DO1575 END DO1576 END DO1577 WHERE (mbathy==0)1578 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1579 END WHERE1580 IF( lk_mpp ) THEN1581 zbathy(:,:) = FLOAT( misfdep(:,:) )1582 CALL lbc_lnk( zbathy, 'T', 1. )1583 misfdep(:,:) = INT( zbathy(:,:) )1584 CALL lbc_lnk( risfdep, 'T', 1. )1585 CALL lbc_lnk( bathy, 'T', 1. )1586 zbathy(:,:) = FLOAT( mbathy(:,:) )1587 CALL lbc_lnk( zbathy, 'T', 1. )1588 mbathy(:,:) = INT( zbathy(:,:) )1589 ENDIF1590 1591 ! fill hole in ice shelf1592 zmisfdep = misfdep1593 zrisfdep = risfdep1594 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1595 DO jj = 2, jpjm11596 DO ji = 2, jpim11597 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1598 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1599 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1600 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1601 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1602 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1603 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1604 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1605 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1606 END IF1607 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1608 misfdep(ji,jj) = ibtest1609 risfdep(ji,jj) = gdepw_1d(ibtest)1610 ENDIF1611 ENDDO1612 ENDDO1613 1614 IF( lk_mpp ) THEN1615 zbathy(:,:) = FLOAT( misfdep(:,:) )1616 CALL lbc_lnk( zbathy, 'T', 1. )1617 misfdep(:,:) = INT( zbathy(:,:) )1618 CALL lbc_lnk( risfdep, 'T', 1. )1619 CALL lbc_lnk( bathy, 'T', 1. )1620 zbathy(:,:) = FLOAT( mbathy(:,:) )1621 CALL lbc_lnk( zbathy, 'T', 1. )1622 mbathy(:,:) = INT( zbathy(:,:) )1623 ENDIF1624 !1625 !! fill hole in bathymetry1626 zmbathy (:,:)=mbathy (:,:)1627 DO jj = 2, jpjm11628 DO ji = 2, jpim11629 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1630 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1631 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1632 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01633 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01634 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01635 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1636 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1637 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1638 END IF1639 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1640 mbathy(ji,jj) = ibtest1641 bathy(ji,jj) = gdepw_1d(ibtest+1)1642 ENDIF1643 END DO1644 END DO1645 IF( lk_mpp ) THEN1646 zbathy(:,:) = FLOAT( misfdep(:,:) )1647 CALL lbc_lnk( zbathy, 'T', 1. )1648 misfdep(:,:) = INT( zbathy(:,:) )1649 CALL lbc_lnk( risfdep, 'T', 1. )1650 CALL lbc_lnk( bathy, 'T', 1. )1651 zbathy(:,:) = FLOAT( mbathy(:,:) )1652 CALL lbc_lnk( zbathy, 'T', 1. )1653 mbathy(:,:) = INT( zbathy(:,:) )1654 ENDIF1655 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1656 DO jj = 1, jpjm11657 DO ji = 1, jpim11658 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1659 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1660 END IF1661 END DO1662 END DO1663 IF( lk_mpp ) THEN1664 zbathy(:,:) = FLOAT( misfdep(:,:) )1665 CALL lbc_lnk( zbathy, 'T', 1. )1666 misfdep(:,:) = INT( zbathy(:,:) )1667 CALL lbc_lnk( risfdep, 'T', 1. )1668 CALL lbc_lnk( bathy, 'T', 1. )1669 zbathy(:,:) = FLOAT( mbathy(:,:) )1670 CALL lbc_lnk( zbathy, 'T', 1. )1671 mbathy(:,:) = INT( zbathy(:,:) )1672 ENDIF1673 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1674 DO jj = 1, jpjm11675 DO ji = 1, jpim11676 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1677 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1678 END IF1679 END DO1680 END DO1681 IF( lk_mpp ) THEN1682 zbathy(:,:) = FLOAT( misfdep(:,:) )1683 CALL lbc_lnk( zbathy, 'T', 1. )1684 misfdep(:,:) = INT( zbathy(:,:) )1685 CALL lbc_lnk( risfdep, 'T', 1. )1686 CALL lbc_lnk( bathy, 'T', 1. )1687 zbathy(:,:) = FLOAT( mbathy(:,:) )1688 CALL lbc_lnk( zbathy, 'T', 1. )1689 mbathy(:,:) = INT( zbathy(:,:) )1690 ENDIF1691 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1692 DO jj = 1, jpjm11693 DO ji = 1, jpi1694 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1695 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1696 END IF1697 END DO1698 END DO1699 IF( lk_mpp ) THEN1700 zbathy(:,:) = FLOAT( misfdep(:,:) )1701 CALL lbc_lnk( zbathy, 'T', 1. )1702 misfdep(:,:) = INT( zbathy(:,:) )1703 CALL lbc_lnk( risfdep, 'T', 1. )1704 CALL lbc_lnk( bathy, 'T', 1. )1705 zbathy(:,:) = FLOAT( mbathy(:,:) )1706 CALL lbc_lnk( zbathy, 'T', 1. )1707 mbathy(:,:) = INT( zbathy(:,:) )1708 ENDIF1709 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1710 DO jj = 1, jpjm11711 DO ji = 1, jpi1712 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1713 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1714 END IF1715 END DO1716 END DO1717 IF( lk_mpp ) THEN1718 zbathy(:,:) = FLOAT( misfdep(:,:) )1719 CALL lbc_lnk( zbathy, 'T', 1. )1720 misfdep(:,:) = INT( zbathy(:,:) )1721 CALL lbc_lnk( risfdep, 'T', 1. )1722 CALL lbc_lnk( bathy, 'T', 1. )1723 zbathy(:,:) = FLOAT( mbathy(:,:) )1724 CALL lbc_lnk( zbathy, 'T', 1. )1725 mbathy(:,:) = INT( zbathy(:,:) )1726 ENDIF1727 ! if not compatible after all check, mask T1728 DO jj = 1, jpj1729 DO ji = 1, jpi1730 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1731 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1732 END IF1733 END DO1734 END DO1735 1736 WHERE (mbathy(:,:) == 1)1737 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1738 END WHERE1739 END DO1740 ! end check compatibility ice shelf/bathy1741 ! remove very shallow ice shelf (less than ~ 10m if 75L)1742 WHERE (misfdep(:,:) <= 5)1743 misfdep = 1; risfdep = 0.0_wp;1744 END WHERE1745 1746 IF( icompt == 0 ) THEN1747 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1748 ELSE1749 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1750 ENDIF1751 1829 1752 1830 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1753 1831 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1754 1755 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1756 1757 END SUBROUTINE 1832 ! 1833 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1834 ! 1835 END SUBROUTINE zgr_isf 1836 1758 1837 1759 1838 SUBROUTINE zgr_sco … … 1801 1880 !! 1802 1881 !!---------------------------------------------------------------------- 1803 !1804 1882 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1805 1883 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1810 1888 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1811 1889 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1812 1890 !! 1813 1891 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1814 1892 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1815 1893 !!---------------------------------------------------------------------- 1816 1894 ! 1817 1895 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1818 1896 ! 1819 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1897 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1820 1898 ! 1821 1899 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1876 1954 DO jj = 1, jpj 1877 1955 DO ji = 1, jpi 1878 IF( bathy(ji,jj) == 0._wp ) THEN 1879 iip1 = MIN( ji+1, jpi ) 1880 ijp1 = MIN( jj+1, jpj ) 1881 iim1 = MAX( ji-1, 1 ) 1882 ijm1 = MAX( jj-1, 1 ) 1883 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1884 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1885 zenv(ji,jj) = rn_sbot_min 1886 ENDIF 1956 IF( bathy(ji,jj) == 0._wp ) THEN 1957 iip1 = MIN( ji+1, jpi ) 1958 ijp1 = MIN( jj+1, jpj ) 1959 iim1 = MAX( ji-1, 1 ) 1960 ijm1 = MAX( jj-1, 1 ) 1961 !!gm BUG fix see ticket #1617 1962 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1963 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1964 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) zenv(ji,jj) = rn_sbot_min 1965 !!gm 1966 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1967 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1968 !!gm zenv(ji,jj) = rn_sbot_min 1969 !!gm ENDIF 1970 !!gm end 1887 1971 ENDIF 1888 1972 END DO … … 1975 2059 ENDIF 1976 2060 ! 1977 IF(lwp) THEN ! Control print1978 WRITE(numout,*)1979 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1980 WRITE(numout,*)1981 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1982 IF( nprint == 1 ) THEN1983 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1984 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1985 ENDIF1986 ENDIF1987 1988 2061 ! ! ============================== 1989 2062 ! ! hbatu, hbatv, hbatf fields … … 2080 2153 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2081 2154 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2082 2083 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2084 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2085 ! 2086 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2087 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2088 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2089 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2090 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2091 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2092 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2155 ! 2156 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2157 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2158 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2159 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2160 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2161 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2162 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2093 2163 2094 2164 #if defined key_agrif 2095 ! Ensure meaningful vertical scale factors in ghost lines/columns 2096 IF( .NOT. Agrif_Root() ) THEN 2097 ! 2098 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2099 e3u_0(1,:,:) = e3u_0(2,:,:) 2100 ENDIF 2101 ! 2102 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2103 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2104 ENDIF 2105 ! 2106 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2107 e3v_0(:,1,:) = e3v_0(:,2,:) 2108 ENDIF 2109 ! 2110 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2111 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2112 ENDIF 2113 ! 2114 ENDIF 2165 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2166 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2167 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2168 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2169 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2170 ENDIF 2115 2171 #endif 2116 2172 2117 fsdept(:,:,:) = gdept_0 (:,:,:) 2118 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2119 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2120 fse3t (:,:,:) = e3t_0 (:,:,:) 2121 fse3u (:,:,:) = e3u_0 (:,:,:) 2122 fse3v (:,:,:) = e3v_0 (:,:,:) 2123 fse3f (:,:,:) = e3f_0 (:,:,:) 2124 fse3w (:,:,:) = e3w_0 (:,:,:) 2125 fse3uw(:,:,:) = e3uw_0 (:,:,:) 2126 fse3vw(:,:,:) = e3vw_0 (:,:,:) 2173 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2174 !!gm and only that !!!!! 2175 !!gm THIS should be removed from here ! 2176 gdept_n(:,:,:) = gdept_0(:,:,:) 2177 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2178 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2179 e3t_n (:,:,:) = e3t_0 (:,:,:) 2180 e3u_n (:,:,:) = e3u_0 (:,:,:) 2181 e3v_n (:,:,:) = e3v_0 (:,:,:) 2182 e3f_n (:,:,:) = e3f_0 (:,:,:) 2183 e3w_n (:,:,:) = e3w_0 (:,:,:) 2184 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2185 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2186 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2187 !! gm end 2127 2188 !! 2128 2189 ! HYBRID : … … 2130 2191 DO ji = 1, jpi 2131 2192 DO jk = 1, jpkm1 2132 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )2133 END DO 2134 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 02193 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2194 END DO 2195 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2135 2196 END DO 2136 2197 END DO … … 2141 2202 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2142 2203 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2143 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2144 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2145 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2146 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2204 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2205 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2206 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2207 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2147 2208 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2148 2209 2149 2210 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2150 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2151 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2152 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2153 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2211 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2212 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2213 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2214 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2154 2215 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2155 2216 ENDIF … … 2183 2244 END DO 2184 2245 ENDIF 2185 2186 !================================================================================2187 ! check the coordinate makes sense2188 !================================================================================2246 ! 2247 !================================================================================ 2248 ! check the coordinate makes sense 2249 !================================================================================ 2189 2250 DO ji = 1, jpi 2190 2251 DO jj = 1, jpj 2191 2252 ! 2192 2253 IF( hbatt(ji,jj) > 0._wp) THEN 2193 2254 DO jk = 1, mbathy(ji,jj) 2194 2255 ! check coordinate is monotonically increasing 2195 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2256 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2196 2257 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2197 2258 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2198 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2199 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2259 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2260 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2200 2261 CALL ctl_stop( ctmp1 ) 2201 2262 ENDIF 2202 2263 ! and check it has never gone negative 2203 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2264 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2204 2265 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2205 2266 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2206 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2207 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2267 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2268 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2208 2269 CALL ctl_stop( ctmp1 ) 2209 2270 ENDIF 2210 2271 ! and check it never exceeds the total depth 2211 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2272 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2212 2273 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2213 2274 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2214 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2275 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2215 2276 CALL ctl_stop( ctmp1 ) 2216 2277 ENDIF 2217 2278 END DO 2218 2279 ! 2219 2280 DO jk = 1, mbathy(ji,jj)-1 2220 2281 ! and check it never exceeds the total depth 2221 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2282 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2222 2283 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2223 2284 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2224 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2285 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2225 2286 CALL ctl_stop( ctmp1 ) 2226 2287 ENDIF 2227 2288 END DO 2228 2229 2289 ENDIF 2230 2231 2290 END DO 2232 2291 END DO … … 2238 2297 END SUBROUTINE zgr_sco 2239 2298 2240 !!====================================================================== 2299 2241 2300 SUBROUTINE s_sh94() 2242 2243 2301 !!---------------------------------------------------------------------- 2244 2302 !! *** ROUTINE s_sh94 *** … … 2251 2309 !! Reference : Song and Haidvogel 1994. 2252 2310 !!---------------------------------------------------------------------- 2253 !2254 2311 INTEGER :: ji, jj, jk ! dummy loop argument 2255 2312 REAL(wp) :: zcoeft, zcoefw ! temporary scalars … … 2257 2314 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2258 2315 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2259 2260 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2261 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2316 !!---------------------------------------------------------------------- 2317 2318 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2319 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2262 2320 2263 2321 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2265 2323 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2266 2324 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2267 2325 ! 2268 2326 DO ji = 1, jpi 2269 2327 DO jj = 1, jpj 2270 2328 ! 2271 2329 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2272 2330 DO jk = 1, jpk … … 2297 2355 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2298 2356 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2299 gdept_0 2300 gdepw_0 2301 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2357 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2358 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2359 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2302 2360 END DO 2303 2361 ! … … 2331 2389 END DO 2332 2390 END DO 2333 2334 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2335 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2336 2391 ! 2392 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2393 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2394 ! 2337 2395 END SUBROUTINE s_sh94 2338 2396 2397 2339 2398 SUBROUTINE s_sf12 2340 2341 2399 !!---------------------------------------------------------------------- 2342 2400 !! *** ROUTINE s_sf12 *** … … 2354 2412 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2355 2413 !!---------------------------------------------------------------------- 2356 !2357 2414 INTEGER :: ji, jj, jk ! dummy loop argument 2358 2415 REAL(wp) :: zsmth ! smoothing around critical depth … … 2361 2418 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2362 2419 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2363 2420 !!---------------------------------------------------------------------- 2364 2421 ! 2365 2422 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2425 2482 2426 2483 DO jk = 1, jpk 2427 gdept_0 2428 gdepw_0 2429 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2484 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2485 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2486 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2430 2487 END DO 2431 2488 … … 2454 2511 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2455 2512 ! 2456 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2513 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2457 2514 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2458 2515 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2467 2524 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2468 2525 ! 2469 ! ! ============= 2470 2471 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2472 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2473 2526 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2527 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2528 ! 2474 2529 END SUBROUTINE s_sf12 2475 2530 2531 2476 2532 SUBROUTINE s_tanh() 2477 2478 2533 !!---------------------------------------------------------------------- 2479 2534 !! *** ROUTINE s_tanh*** … … 2485 2540 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2486 2541 !!---------------------------------------------------------------------- 2487 2488 INTEGER :: ji, jj, jk ! dummy loop argument 2542 INTEGER :: ji, jj, jk ! dummy loop argument 2489 2543 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2490 2491 2544 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2492 2545 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2493 2494 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2495 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2546 !!---------------------------------------------------------------------- 2547 2548 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2549 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2496 2550 2497 2551 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2523 2577 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2524 2578 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2525 gdept_0 2526 gdepw_0 2527 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2579 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2580 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2581 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2528 2582 END DO 2529 2583 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2542 2596 END DO 2543 2597 END DO 2544 2545 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2546 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2547 2598 ! 2599 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2600 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2601 ! 2548 2602 END SUBROUTINE s_tanh 2603 2549 2604 2550 2605 FUNCTION fssig( pk ) RESULT( pf ) … … 2618 2673 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2619 2674 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2620 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2621 REAL(wp), INTENT(in ) :: pzs ! surface box depth2622 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2623 REAL(wp) :: za1,za2,za3 ! local variables2624 REAL(wp) :: zn1,zn2 ! local variables2625 REAL(wp) :: za,zb,zx ! local variables2626 integer :: jk2627 !!----------------------------------------------------------------------2628 ! 2629 2630 zn1 = 1. /(jpk-1.)2631 zn2 = 1. - zn12632 2675 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2676 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2677 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2678 ! 2679 INTEGER :: jk ! dummy loop index 2680 REAL(wp) :: za1,za2,za3 ! local scalar 2681 REAL(wp) :: zn1,zn2 ! - - 2682 REAL(wp) :: za,zb,zx ! - - 2683 !!---------------------------------------------------------------------- 2684 ! 2685 zn1 = 1._wp / REAL( jpkm1, wp ) 2686 zn2 = 1._wp - zn1 2687 ! 2633 2688 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2634 2689 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2635 2690 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2636 2691 ! 2637 2692 za = pzb - za3*(pzs-za1)-za2 2638 2693 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2639 2694 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2640 2695 zx = 1.0_wp-za/2.0_wp-zb 2641 2696 ! 2642 2697 DO jk = 1, jpk 2643 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2644 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2645 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2698 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2699 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2700 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2646 2701 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2647 ENDDO 2648 2702 END DO 2649 2703 ! 2650 2704 END FUNCTION fgamma
Note: See TracChangeset
for help on using the changeset viewer.
