- Timestamp:
- 2015-10-16T19:19:43+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
r5790 r5802 64 64 REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio 65 65 !! 66 REAL(wp), DIMENSION(:,: ), POINTER :: zde3t 66 REAL(wp), DIMENSION(:,: ), POINTER :: zde3t, zdtem, zdsal 67 67 REAL(wp), DIMENSION(:,: ), POINTER :: zssh0 68 68 REAL(wp), DIMENSION(:,:,: ), POINTER :: ztmp3d … … 74 74 75 75 CALL wrk_alloc(jpi,jpj,jpk, ztmp3d ) 76 CALL wrk_alloc(jpi,jpj, zde3t )76 CALL wrk_alloc(jpi,jpj, zde3t , zdtem, zdsal ) 77 77 CALL wrk_alloc(jpi,jpj, zssh0 ) 78 78 79 79 ! get unbalance (volume heat and salt) 80 80 ! initialisation 81 zde3t (:,:) = 0.0_wp 82 pvol_flx(:,:,: ) = 0.0_wp 83 pts_flx (:,:,:,:) = 0.0_wp 84 85 zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt 86 IF (lwp) PRINT *, 'total volume correction 0 = ',zsum 87 zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 88 IF (lwp) PRINT *, 'total heat correction 0 = ',zsum 89 zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 90 IF (lwp) PRINT *, 'total salt correction 0 = ',zsum 91 92 ! mask tsn and tsb (should be useless) 93 tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); 94 tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); 95 96 ! diagnose non conservation of heat, salt and volume 97 r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 98 zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 99 IF ( lk_vvl ) zssh0 = 0.0_wp 100 DO jk = 1,jpk-1 101 DO ji = 2,jpi-1 102 DO jj = 2,jpj-1 103 ! volume differences 104 zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); 105 106 ! shh changes 107 IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 108 zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 109 zssh0(ji,jj) = 0._wp 110 END IF 111 112 ! ocean cell now 113 ! case where we open, enlarge or thin a cell : 114 pvol_flx(ji,jj,jk) = zde3t(ji,jj) * r1_tiscpl 115 pts_flx (ji,jj,jk,jp_sal)= tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 116 pts_flx (ji,jj,jk,jp_tem)= tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl 117 END DO 118 END DO 119 END DO 120 ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0 121 PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt 122 zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt 123 IF (lwp) PRINT *, 'total volume correction 1 = ',zsum 124 zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 125 IF (lwp) PRINT *, 'total heat correction 1 = ',zsum 126 zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 127 IF (lwp) PRINT *, 'total salt correction 1 = ',zsum 128 129 zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 130 IF ( lk_vvl ) zssh0 = 0.0_wp 131 DO jk = 1,jpk-1 132 DO ji = 2,jpi-1 133 DO jj = 2,jpj-1 134 ! volume differences 135 zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); 136 137 ! shh changes 138 IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 139 zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 140 zssh0(ji,jj) = 0._wp 141 END IF 142 143 ! ocean cell before and mask cell now 144 IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN 145 ! case where we close a cell and adjacent cell open 146 pvol_flx(ji,jj,jk) = zde3t(ji,jj) * r1_tiscpl 147 pts_flx (ji,jj,jk,jp_sal)= tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 148 pts_flx (ji,jj,jk,jp_tem)= tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl 149 150 jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ; 151 152 zsum = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) + e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) & 153 & + e12t(jim1,jj ) * tmask(jim1,jj ,jk) + e12t(jip1,jj ) * tmask(jip1,jj ,jk) 154 155 IF ( zsum .NE. 0._wp ) THEN 156 zjip1_ratio = e12t(jip1,jj ) * tmask(jip1,jj ,jk) / zsum 157 zjim1_ratio = e12t(jim1,jj ) * tmask(jim1,jj ,jk) / zsum 158 zjjp1_ratio = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) / zsum 159 zjjm1_ratio = e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) / zsum 160 161 pvol_flx(ji ,jjp1,jk ) = pvol_flx(ji ,jjp1,jk ) + pvol_flx(ji,jj,jk ) * zjjp1_ratio 162 pvol_flx(ji ,jjm1,jk ) = pvol_flx(ji ,jjm1,jk ) + pvol_flx(ji,jj,jk ) * zjjm1_ratio 163 pvol_flx(jip1,jj ,jk ) = pvol_flx(jip1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjip1_ratio 164 pvol_flx(jim1,jj ,jk ) = pvol_flx(jim1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjim1_ratio 165 pts_flx (ji ,jjp1,jk,jp_sal) = pts_flx (ji ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio 166 pts_flx (ji ,jjm1,jk,jp_sal) = pts_flx (ji ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio 167 pts_flx (jip1,jj ,jk,jp_sal) = pts_flx (jip1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio 168 pts_flx (jim1,jj ,jk,jp_sal) = pts_flx (jim1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio 169 pts_flx (ji ,jjp1,jk,jp_tem) = pts_flx (ji ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio 170 pts_flx (ji ,jjm1,jk,jp_tem) = pts_flx (ji ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio 171 pts_flx (jip1,jj ,jk,jp_tem) = pts_flx (jip1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio 172 pts_flx (jim1,jj ,jk,jp_tem) = pts_flx (jim1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio 173 174 ! set to 0 the cell we distributed over neigbourg cells 175 pvol_flx(ji,jj,jk ) = 0._wp 176 pts_flx (ji,jj,jk,jp_sal) = 0._wp 177 pts_flx (ji,jj,jk,jp_tem) = 0._wp 178 179 ELSE IF (zsum .EQ. 0._wp ) THEN 180 ! case where we close a cell and no adjacent cell open 181 ! check if the cell beneath is wet 182 IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN 183 pvol_flx(ji,jj,jk+1) = pvol_flx(ji,jj,jk+1) + pvol_flx(ji,jj,jk) 184 pts_flx (ji,jj,jk+1,jp_sal)= pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal) 185 pts_flx (ji,jj,jk+1,jp_tem)= pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem) 186 187 ! set to 0 the cell we distributed over neigbourg cells 188 pvol_flx(ji,jj,jk ) = 0._wp 189 pts_flx (ji,jj,jk,jp_sal) = 0._wp 190 pts_flx (ji,jj,jk,jp_tem) = 0._wp 191 ELSE 192 ! case no adjacent cell on the horizontal and on the vertical 193 PRINT *, 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' 194 PRINT *, ' ',mig(ji),' ',mjg(jj),' ',jk 195 PRINT *, ' ',ji,' ',jj,' ',jk,' ',narea 196 PRINT *, ' we are now looking for the closest wet cell on the horizontal ' 197 ! We deal with this points later. 198 END IF 199 END IF 200 END IF 201 END DO 202 END DO 203 END DO 204 205 zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt 206 IF (lwp) PRINT *, 'total volume correction 2 = ',zsum 207 zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 208 IF (lwp) PRINT *, 'total heat correction 2 = ',zsum 209 zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 210 IF (lwp) PRINT *, 'total salt correction 2 = ',zsum 211 212 ! allocation and initialisation of the list of problematic point 213 ALLOCATE(vnpts(jpnij)) 214 vnpts(:)=0 215 216 ! fill narea location with the number of problematic point 217 DO jk = 1,jpk-1 218 DO ji = 2,jpi-1 219 DO jj = 2,jpj-1 220 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 221 & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 ) THEN 222 vnpts(narea) = vnpts(narea) + 1 223 END IF 224 END DO 225 END DO 226 END DO 227 228 ! build array of total problematic point on each cpu (share to each cpu) 229 CALL mpp_max(vnpts,jpnij) 230 231 ! size of the new variable 232 npts = SUM(vnpts) 233 234 ! allocation of the coordinates, correction, index vector for the problematic points 235 ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) 236 ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20 237 zcorr_vol(:) = 0.0_wp 238 zcorr_sal(:) = 0.0_wp 239 zcorr_tem(:) = 0.0_wp 240 241 ! fill new variable 242 jpts = SUM(vnpts(1:narea-1)) 243 DO jk = 1,jpk-1 244 DO ji = 2,jpi-1 245 DO jj = 2,jpj-1 246 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 247 & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 ) THEN 248 jpts = jpts + 1 ! positioning in the vnpts vector for the area narea 249 PRINT *, 'corrected point ', narea, ji, jj, jk, jpts 250 ixpts(jpts) = ji ; iypts(jpts) = jj ; izpts(jpts) = jk 251 zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) 252 zcorr_vol(jpts) = pvol_flx(ji,jj,jk) 253 zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) 254 zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) 255 ! set flx to 0 (safer) 256 pvol_flx(ji,jj,jk ) = 0.0_wp 257 pts_flx (ji,jj,jk,jp_sal) = 0.0_wp 258 pts_flx (ji,jj,jk,jp_tem) = 0.0_wp 259 PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt 260 END IF 261 END DO 262 END DO 263 END DO 264 265 ! build array of total problematic point on each cpu (share to each cpu) 266 CALL mpp_max(zlat ,npts) 267 CALL mpp_max(zlon ,npts) 268 CALL mpp_max(izpts,npts) 269 270 ! put correction term in the closest cell 271 PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts 272 DO jpts = 1,npts 273 CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts)) 274 PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts) 275 DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) 276 DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) 277 jk = izpts(jpts) 278 pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk ) + zcorr_vol(jpts) 279 pts_flx (ji,jj,jk,jp_sal) = pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) 280 pts_flx (ji,jj,jk,jp_tem) = pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) 281 END DO 282 END DO 283 END DO 284 ! deallocate variables 285 DEALLOCATE(vnpts) 286 DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat) 287 288 ! add contribution store on the hallo (lbclnk remove one of the contribution) 289 pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:) 290 pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) 291 pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:) 292 293 CALL lbc_sum(pvol_flx(:,:,: ),'T',1.) 294 CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 295 CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 296 297 ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!! 298 zsumn = glob_sum ( fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt 299 ztmp3d(:,:,:) = 0.0 300 ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 301 zsumb = glob_sum_full(ztmp3d) 302 zsum = glob_sum ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt) 303 IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum 304 ! CHECK salt 305 zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 306 zsumb = glob_sum( tsb(:,:,:,jp_sal) * pe3t_b(:,:,:) * ptmask_b(:,:,:)) 307 zsum = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt) 308 IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum 309 ! CHECK heat 310 zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 311 zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * ptmask_b(:,:,:)) 312 zsum = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt) 313 IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum 314 !! 315 CALL wrk_dealloc(jpi,jpj,jpk, ztmp3d ) 316 CALL wrk_dealloc(jpi,jpj, zde3t ) 317 CALL wrk_dealloc(jpi,jpj, zssh0 ) 81 zde3t (:,:) = 0.0_wp 82 pvol_flx(:,:,: ) = 0.0_wp 83 pts_flx (:,:,:,:) = 0.0_wp 84 85 ! mask tsn and tsb (should be useless) 86 tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); 87 tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); 88 89 ! diagnose non conservation of heat, salt and volume 90 r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 91 92 zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 93 IF ( lk_vvl ) zssh0 = 0.0_wp ! already include in the levels by definition 94 95 DO jk = 1,jpk-1 96 DO ji = 2,jpi-1 97 DO jj = 2,jpj-1 98 IF (tmask_h(ji,jj) == 1._wp) THEN 99 100 ! volume differences 101 zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk) 102 103 ! heat diff 104 zdtem(ji,jj) = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) & 105 - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 106 ! salt diff 107 zdsal(ji,jj) = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) & 108 - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 109 110 ! shh changes 111 IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 112 zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 113 zssh0(ji,jj) = 0._wp 114 END IF 115 116 ! volume, heat and salt differences in each cell 117 pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk) + zde3t(ji,jj) * r1_tiscpl 118 pts_flx (ji,jj,jk,jp_sal)= pts_flx (ji,jj,jk,jp_sal) + zdsal(ji,jj) * r1_tiscpl 119 pts_flx (ji,jj,jk,jp_tem)= pts_flx (ji,jj,jk,jp_tem) + zdtem(ji,jj) * r1_tiscpl 120 121 IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN 122 ! case where we close a cell: check if the neighbour cells are wet 123 124 jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ; 125 126 zsum = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) + e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) & 127 & + e12t(jim1,jj ) * tmask(jim1,jj ,jk) + e12t(jip1,jj ) * tmask(jip1,jj ,jk) 128 129 IF ( zsum .NE. 0._wp ) THEN 130 zjip1_ratio = e12t(jip1,jj ) * tmask(jip1,jj ,jk) / zsum 131 zjim1_ratio = e12t(jim1,jj ) * tmask(jim1,jj ,jk) / zsum 132 zjjp1_ratio = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) / zsum 133 zjjm1_ratio = e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) / zsum 134 135 pvol_flx(ji ,jjp1,jk ) = pvol_flx(ji ,jjp1,jk ) + pvol_flx(ji,jj,jk ) * zjjp1_ratio 136 pvol_flx(ji ,jjm1,jk ) = pvol_flx(ji ,jjm1,jk ) + pvol_flx(ji,jj,jk ) * zjjm1_ratio 137 pvol_flx(jip1,jj ,jk ) = pvol_flx(jip1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjip1_ratio 138 pvol_flx(jim1,jj ,jk ) = pvol_flx(jim1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjim1_ratio 139 pts_flx (ji ,jjp1,jk,jp_sal) = pts_flx (ji ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio 140 pts_flx (ji ,jjm1,jk,jp_sal) = pts_flx (ji ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio 141 pts_flx (jip1,jj ,jk,jp_sal) = pts_flx (jip1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio 142 pts_flx (jim1,jj ,jk,jp_sal) = pts_flx (jim1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio 143 pts_flx (ji ,jjp1,jk,jp_tem) = pts_flx (ji ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio 144 pts_flx (ji ,jjm1,jk,jp_tem) = pts_flx (ji ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio 145 pts_flx (jip1,jj ,jk,jp_tem) = pts_flx (jip1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio 146 pts_flx (jim1,jj ,jk,jp_tem) = pts_flx (jim1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio 147 148 ! set to 0 the cell we distributed over neigbourg cells 149 pvol_flx(ji,jj,jk ) = 0._wp 150 pts_flx (ji,jj,jk,jp_sal) = 0._wp 151 pts_flx (ji,jj,jk,jp_tem) = 0._wp 152 153 ELSE IF (zsum .EQ. 0._wp ) THEN 154 ! case where we close a cell and no adjacent cell open 155 ! check if the cell beneath is wet 156 IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN 157 pvol_flx(ji,jj,jk+1) = pvol_flx(ji,jj,jk+1) + pvol_flx(ji,jj,jk) 158 pts_flx (ji,jj,jk+1,jp_sal)= pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal) 159 pts_flx (ji,jj,jk+1,jp_tem)= pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem) 160 161 ! set to 0 the cell we distributed over neigbourg cells 162 pvol_flx(ji,jj,jk ) = 0._wp 163 pts_flx (ji,jj,jk,jp_sal) = 0._wp 164 pts_flx (ji,jj,jk,jp_tem) = 0._wp 165 ELSE 166 ! case no adjacent cell on the horizontal and on the vertical 167 PRINT *, 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' 168 PRINT *, ' ',mig(ji),' ',mjg(jj),' ',jk 169 PRINT *, ' ',ji,' ',jj,' ',jk,' ',narea 170 PRINT *, ' we are now looking for the closest wet cell on the horizontal ' 171 ! We deal with this points later. 172 END IF 173 END IF 174 END IF 175 END IF 176 END DO 177 END DO 178 END DO 179 180 CALL lbc_sum(pvol_flx(:,:,: ),'T',1.) 181 CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 182 CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 183 184 zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt 185 IF (lwp) PRINT *, 'total volume correction 21 = ',zsum 186 zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 187 IF (lwp) PRINT *, 'total heat correction 21 = ',zsum 188 zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 189 IF (lwp) PRINT *, 'total salt correction 21 = ',zsum 190 191 ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point 192 ! allocation and initialisation of the list of problematic point 193 ALLOCATE(vnpts(jpnij)) 194 vnpts(:)=0 195 196 ! fill narea location with the number of problematic point 197 DO jk = 1,jpk-1 198 DO ji = 2,jpi-1 199 DO jj = 2,jpj-1 200 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 201 & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 & 202 & .AND. tmask_h(ji,jj) == 1._wp ) THEN 203 vnpts(narea) = vnpts(narea) + 1 204 END IF 205 END DO 206 END DO 207 END DO 208 209 ! build array of total problematic point on each cpu (share to each cpu) 210 CALL mpp_max(vnpts,jpnij) 211 212 ! size of the new variable 213 npts = SUM(vnpts) 214 215 ! allocation of the coordinates, correction, index vector for the problematic points 216 ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) 217 ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20 218 zcorr_vol(:) = 0.0_wp 219 zcorr_sal(:) = 0.0_wp 220 zcorr_tem(:) = 0.0_wp 221 222 ! fill new variable 223 jpts = SUM(vnpts(1:narea-1)) 224 DO jk = 1,jpk-1 225 DO ji = 2,jpi-1 226 DO jj = 2,jpj-1 227 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 228 & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 & 229 & .AND. tmask_h(ji,jj) == 1 ) THEN 230 jpts = jpts + 1 ! positioning in the vnpts vector for the area narea 231 PRINT *, 'corrected point ', narea, ji, jj, jk, jpts 232 ixpts(jpts) = ji ; iypts(jpts) = jj ; izpts(jpts) = jk 233 zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) 234 zcorr_vol(jpts) = pvol_flx(ji,jj,jk) 235 zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) 236 zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) 237 ! set flx to 0 (safer) 238 pvol_flx(ji,jj,jk ) = 0.0_wp 239 pts_flx (ji,jj,jk,jp_sal) = 0.0_wp 240 pts_flx (ji,jj,jk,jp_tem) = 0.0_wp 241 PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt 242 END IF 243 END DO 244 END DO 245 END DO 246 247 ! build array of total problematic point on each cpu (share to each cpu) 248 CALL mpp_max(zlat ,npts) 249 CALL mpp_max(zlon ,npts) 250 CALL mpp_max(izpts,npts) 251 252 ! put correction term in the closest cell 253 PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts 254 DO jpts = 1,npts 255 CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts)) 256 PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts) 257 DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) 258 DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) 259 jk = izpts(jpts) 260 pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk ) + zcorr_vol(jpts) 261 pts_flx (ji,jj,jk,jp_sal) = pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) 262 pts_flx (ji,jj,jk,jp_tem) = pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) 263 END DO 264 END DO 265 END DO 266 ! deallocate variables 267 DEALLOCATE(vnpts) 268 DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat) 269 270 ! add contribution store on the hallo (lbclnk remove one of the contribution) 271 pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:) 272 pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:) 273 pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) 274 275 CALL lbc_sum(pvol_flx(:,:,: ),'T',1.) 276 CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 277 CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 278 279 ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!! 280 zsum = glob_sum_full( pvol_flx(:,:,:) ) 281 IF (lwp) PRINT *, 'CHECK vol = ',zsum 282 ! CHECK salt 283 zsum = glob_sum( pts_flx(:,:,:,jp_sal) ) 284 IF (lwp) PRINT *, 'CHECK salt = ',zsum 285 ! CHECK heat 286 zsum = glob_sum( pts_flx(:,:,:,jp_tem) ) 287 IF (lwp) PRINT *, 'CHECK heat = ',zsum 288 !! 289 CALL wrk_dealloc(jpi,jpj,jpk, ztmp3d ) 290 CALL wrk_dealloc(jpi,jpj, zde3t ) 291 CALL wrk_dealloc(jpi,jpj, zssh0 ) 318 292 END SUBROUTINE iscpl_cons 319 293
Note: See TracChangeset
for help on using the changeset viewer.
