Changeset 11771
- Timestamp:
- 2019-10-23T15:04:25+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC
- Files:
-
- 1 added
- 5 edited
-
diawri.F90 (modified) (48 diffs)
-
domvvl.F90 (modified) (35 diffs)
-
sbcmod.F90 (added)
-
stpctl.F90 (modified) (5 diffs)
-
trazdf.F90 (modified) (13 diffs)
-
usrdef_sbc.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90
r10425 r11771 57 57 USE lib_mpp ! MPP library 58 58 USE timing ! preformance summary 59 USE diu rnal_bulk! diurnal warm layer60 USE cool_skin! Cool skin59 USE diu_bulk ! diurnal warm layer 60 USE diu_coolskin ! Cool skin 61 61 62 62 IMPLICIT NONE … … 97 97 98 98 99 SUBROUTINE dia_wri( kt )99 SUBROUTINE dia_wri( kt, Kmm ) 100 100 !!--------------------------------------------------------------------- 101 101 !! *** ROUTINE dia_wri *** … … 107 107 !!---------------------------------------------------------------------- 108 108 INTEGER, INTENT( in ) :: kt ! ocean time-step index 109 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 109 110 !! 110 111 INTEGER :: ji, jj, jk ! dummy loop indices … … 115 116 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 116 117 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace 117 REAL(wp), DIMENSION(jpi,jpj,jpk) :: bu, bv ! volume of u- and v-boxes118 REAL(wp), DIMENSION(jpi,jpj,jpk) :: r1_bt ! inverse of t-box volume119 118 !!---------------------------------------------------------------------- 120 119 ! … … 123 122 ! Output the initial state and forcings 124 123 IF( ninist == 1 ) THEN 125 CALL dia_wri_state( 'output.init' )124 CALL dia_wri_state( Kmm, 'output.init' ) 126 125 ninist = 0 127 126 ENDIF … … 132 131 CALL iom_put("e3v_0", e3v_0(:,:,:) ) 133 132 ! 134 CALL iom_put( "e3t" , e3t _n(:,:,:) )135 CALL iom_put( "e3u" , e3u _n(:,:,:) )136 CALL iom_put( "e3v" , e3v _n(:,:,:) )137 CALL iom_put( "e3w" , e3w _n(:,:,:) )133 CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 134 CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 135 CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 136 CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 138 137 IF( iom_use("e3tdef") ) & 139 CALL iom_put( "e3tdef" , ( ( e3t _n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )138 CALL iom_put( "e3tdef" , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 140 139 141 140 IF( ll_wd ) THEN 142 CALL iom_put( "ssh" , (ssh n+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)141 CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying) 143 142 ELSE 144 CALL iom_put( "ssh" , ssh n) ! sea surface height143 CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height 145 144 ENDIF 146 145 147 146 IF( iom_use("wetdep") ) & ! wet depth 148 CALL iom_put( "wetdep" , ht_0(:,:) + ssh n(:,:) )147 CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 149 148 150 CALL iom_put( "toce", ts n(:,:,:,jp_tem) ) ! 3D temperature151 CALL iom_put( "sst", ts n(:,:,1,jp_tem) ) ! surface temperature149 CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature 150 CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature 152 151 IF ( iom_use("sbt") ) THEN 153 152 DO jj = 1, jpj 154 153 DO ji = 1, jpi 155 154 ikbot = mbkt(ji,jj) 156 z2d(ji,jj) = ts n(ji,jj,ikbot,jp_tem)155 z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 157 156 END DO 158 157 END DO … … 160 159 ENDIF 161 160 162 CALL iom_put( "soce", ts n(:,:,:,jp_sal) ) ! 3D salinity163 CALL iom_put( "sss", ts n(:,:,1,jp_sal) ) ! surface salinity161 CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity 162 CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity 164 163 IF ( iom_use("sbs") ) THEN 165 164 DO jj = 1, jpj 166 165 DO ji = 1, jpi 167 166 ikbot = mbkt(ji,jj) 168 z2d(ji,jj) = ts n(ji,jj,ikbot,jp_sal)167 z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 169 168 END DO 170 169 END DO … … 177 176 DO jj = 2, jpjm1 178 177 DO ji = fs_2, fs_jpim1 ! vector opt. 179 zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * u n(ji ,jj,mbku(ji ,jj)) )**2 &180 & + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * u n(ji-1,jj,mbku(ji-1,jj)) )**2 &181 & + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * v n(ji,jj ,mbkv(ji,jj )) )**2 &182 & + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * v n(ji,jj-1,mbkv(ji,jj-1)) )**2178 zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 & 179 & + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 & 180 & + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 & 181 & + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2 183 182 z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 184 183 ! … … 189 188 ENDIF 190 189 191 CALL iom_put( "uoce", u n(:,:,:) ) ! 3D i-current192 CALL iom_put( "ssu", u n(:,:,1) ) ! surface i-current190 CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current 191 CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current 193 192 IF ( iom_use("sbu") ) THEN 194 193 DO jj = 1, jpj 195 194 DO ji = 1, jpi 196 195 ikbot = mbku(ji,jj) 197 z2d(ji,jj) = u n(ji,jj,ikbot)196 z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 198 197 END DO 199 198 END DO … … 201 200 ENDIF 202 201 203 CALL iom_put( "voce", v n(:,:,:) ) ! 3D j-current204 CALL iom_put( "ssv", v n(:,:,1) ) ! surface j-current202 CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current 203 CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current 205 204 IF ( iom_use("sbv") ) THEN 206 205 DO jj = 1, jpj 207 206 DO ji = 1, jpi 208 207 ikbot = mbkv(ji,jj) 209 z2d(ji,jj) = v n(ji,jj,ikbot)208 z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 210 209 END DO 211 210 END DO … … 213 212 ENDIF 214 213 215 CALL iom_put( "woce", w n) ! vertical velocity214 CALL iom_put( "woce", ww ) ! vertical velocity 216 215 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value 217 216 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 218 217 z2d(:,:) = rau0 * e1e2t(:,:) 219 218 DO jk = 1, jpk 220 z3d(:,:,jk) = w n(:,:,jk) * z2d(:,:)219 z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:) 221 220 END DO 222 221 CALL iom_put( "w_masstr" , z3d ) … … 236 235 DO jj = 2, jpjm1 ! sal gradient 237 236 DO ji = fs_2, fs_jpim1 ! vector opt. 238 zztmp = ts n(ji,jj,jk,jp_sal)239 zztmpx = ( ts n(ji+1,jj,jk,jp_sal) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,jk,jp_sal) ) * r1_e1u(ji-1,jj)240 zztmpy = ( ts n(ji,jj+1,jk,jp_sal) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,jk,jp_sal) ) * r1_e2v(ji,jj-1)237 zztmp = ts(ji,jj,jk,jp_sal,Kmm) 238 zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj) 239 zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1) 241 240 z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) & 242 241 & * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) … … 253 252 DO jj = 2, jpjm1 ! sst gradient 254 253 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zztmp = ts n(ji,jj,1,jp_tem)256 zztmpx = ( ts n(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) * r1_e1u(ji-1,jj)257 zztmpy = ( ts n(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)254 zztmp = ts(ji,jj,1,jp_tem,Kmm) 255 zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 256 zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 258 257 z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) & 259 258 & * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) … … 272 271 DO jj = 1, jpj 273 272 DO ji = 1, jpi 274 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)273 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 275 274 END DO 276 275 END DO … … 284 283 DO jj = 1, jpj 285 284 DO ji = 1, jpi 286 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)285 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 287 286 END DO 288 287 END DO … … 296 295 DO jj = 1, jpj 297 296 DO ji = 1, jpi 298 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)297 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 299 298 END DO 300 299 END DO … … 306 305 z3d(:,:,jpk) = 0._wp 307 306 DO jk = 1, jpkm1 308 DO jj = 2, jpj 309 DO ji = 2, jpi 310 zztmpx = 0.5 * ( un(ji-1,jj ,jk) + un(ji,jj,jk) ) 311 zztmpy = 0.5 * ( vn(ji ,jj-1,jk) + vn(ji,jj,jk) ) 312 z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy ) 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 ! vector opt. 309 zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 310 z3d(ji,jj,jk) = zztmp * ( uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 311 & + uu(ji ,jj,jk,Kmm)**2 * e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 312 & + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 313 & + vv(ji,jj ,jk,Kmm)**2 * e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) ) 313 314 END DO 314 315 END DO … … 326 327 DO jj = 2, jpj 327 328 DO ji = 2, jpi 328 z3d(ji,jj,jk) = 0.25_wp * ( u n(ji ,jj,jk) * un(ji ,jj,jk) * e1e2u(ji ,jj) * e3u_n(ji ,jj,jk) &329 & + u n(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) &330 & + v n(ji,jj ,jk) * vn(ji,jj ,jk) * e1e2v(ji,jj ) * e3v_n(ji,jj ,jk) &331 & + v n(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk) ) &332 & * r1_e1e2t(ji,jj) / e3t _n(ji,jj,jk) * tmask(ji,jj,jk)329 z3d(ji,jj,jk) = 0.25_wp * ( uu(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 330 & + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 331 & + vv(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,jk,Kmm) & 332 & + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) ) & 333 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 333 334 END DO 334 335 END DO … … 342 343 DO jj = 1, jpj 343 344 DO ji = 1, jpi 344 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * z3d(ji,jj,jk) * tmask(ji,jj,jk)345 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 345 346 END DO 346 347 END DO … … 350 351 ENDIF 351 352 ! 352 CALL iom_put( "hdiv", hdiv n) ! Horizontal divergence353 CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence 353 354 354 355 IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN … … 358 359 DO jj = 1, jpjm1 359 360 DO ji = 1, fs_jpim1 ! vector opt. 360 z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &361 & - e1u(ji ,jj+1) * u n(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj)361 z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm) & 362 & - e1u(ji ,jj+1) * uu(ji ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm) ) * r1_e1e2f(ji,jj) 362 363 END DO 363 364 END DO … … 378 379 DO jj = 1, jpjm1 379 380 DO ji = 1, fs_jpim1 ! vector opt. 380 ze3 = ( e3t _n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &381 & + e3t _n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) )381 ze3 = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 382 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 382 383 IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3 383 384 ELSE ; ze3 = 0._wp … … 397 398 z2d(:,:) = 0.e0 398 399 DO jk = 1, jpkm1 399 z3d(:,:,jk) = rau0 * u n(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)400 z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk) 400 401 z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 401 402 END DO … … 409 410 DO jj = 2, jpjm1 410 411 DO ji = fs_2, fs_jpim1 ! vector opt. 411 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts n(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )412 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 412 413 END DO 413 414 END DO … … 422 423 DO jj = 2, jpjm1 423 424 DO ji = fs_2, fs_jpim1 ! vector opt. 424 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts n(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )425 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 425 426 END DO 426 427 END DO … … 434 435 z3d(:,:,jpk) = 0.e0 435 436 DO jk = 1, jpkm1 436 z3d(:,:,jk) = rau0 * v n(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)437 z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk) 437 438 END DO 438 439 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction … … 444 445 DO jj = 2, jpjm1 445 446 DO ji = fs_2, fs_jpim1 ! vector opt. 446 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts n(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )447 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 447 448 END DO 448 449 END DO … … 457 458 DO jj = 2, jpjm1 458 459 DO ji = fs_2, fs_jpim1 ! vector opt. 459 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts n(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )460 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 460 461 END DO 461 462 END DO … … 470 471 DO jj = 2, jpjm1 471 472 DO ji = fs_2, fs_jpim1 ! vector opt. 472 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)473 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 473 474 END DO 474 475 END DO … … 482 483 DO jj = 2, jpjm1 483 484 DO ji = fs_2, fs_jpim1 ! vector opt. 484 z2d(ji,jj) = z2d(ji,jj) + e3t _n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)485 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 485 486 END DO 486 487 END DO … … 493 494 ! 494 495 495 IF (ln_diatmb) CALL dia_tmb 496 IF (ln_diatmb) CALL dia_tmb( Kmm ) ! tmb values 496 497 497 IF (ln_dia25h) CALL dia_25h( kt )! 25h averaging498 IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging 498 499 499 500 IF( ln_timing ) CALL timing_stop('dia_wri') … … 521 522 522 523 523 SUBROUTINE dia_wri( kt )524 SUBROUTINE dia_wri( kt, Kmm ) 524 525 !!--------------------------------------------------------------------- 525 526 !! *** ROUTINE dia_wri *** … … 534 535 !!---------------------------------------------------------------------- 535 536 INTEGER, INTENT( in ) :: kt ! ocean time-step index 537 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 536 538 ! 537 539 LOGICAL :: ll_print = .FALSE. ! =T print and flush numout … … 551 553 ! 552 554 IF( ninist == 1 ) THEN !== Output the initial state and forcings ==! 553 CALL dia_wri_state( 'output.init' )555 CALL dia_wri_state( Kmm, 'output.init' ) 554 556 ninist = 0 555 557 ENDIF … … 688 690 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 689 691 IF( .NOT.ln_linssh ) THEN 690 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t _n692 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t(:,:,:,Kmm) 691 693 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 692 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t _n694 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t(:,:,:,Kmm) 693 695 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 694 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t _n696 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t(:,:,:,Kmm) 695 697 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 696 698 ENDIF … … 709 711 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 710 712 IF( ln_linssh ) THEN 711 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts n(:,:,1,jp_tem)713 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm) 712 714 & , "KgC/m2/s", & ! sosst_cd 713 715 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 714 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts n(:,:,1,jp_sal)716 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm) 715 717 & , "KgPSU/m2/s",& ! sosss_cd 716 718 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) … … 797 799 798 800 ! !!! nid_U : 3D 799 CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! u n801 CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! uu(:,:,:,Kmm) 800 802 & jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 801 803 IF( ln_wave .AND. ln_sdw) THEN … … 810 812 811 813 ! !!! nid_V : 3D 812 CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! v n814 CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vv(:,:,:,Kmm) 813 815 & jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 814 816 IF( ln_wave .AND. ln_sdw) THEN … … 823 825 824 826 ! !!! nid_W : 3D 825 CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! w n827 CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! ww 826 828 & jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 827 829 CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt … … 861 863 862 864 IF( .NOT.ln_linssh ) THEN 863 CALL histwrite( nid_T, "votemper", it, ts n(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content864 CALL histwrite( nid_T, "vosaline", it, ts n(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content865 CALL histwrite( nid_T, "sosstsst", it, ts n(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content866 CALL histwrite( nid_T, "sosaline", it, ts n(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content865 CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T ) ! heat content 866 CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T ) ! salt content 867 CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT ) ! sea surface heat content 868 CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity content 867 869 ELSE 868 CALL histwrite( nid_T, "votemper", it, ts n(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature869 CALL histwrite( nid_T, "vosaline", it, ts n(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity870 CALL histwrite( nid_T, "sosstsst", it, ts n(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature871 CALL histwrite( nid_T, "sosaline", it, ts n(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity870 CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature 871 CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity 872 CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature 873 CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity 872 874 ENDIF 873 875 IF( .NOT.ln_linssh ) THEN 874 zw3d(:,:,:) = ( ( e3t _n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2875 CALL histwrite( nid_T, "vovvle3t", it, e3t _n (:,:,:) , ndim_T , ndex_T ) ! level thickness876 CALL histwrite( nid_T, "vovvldep", it, gdept _n(:,:,:) , ndim_T , ndex_T ) ! t-point depth876 zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 877 CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T ) ! level thickness 878 CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T ) ! t-point depth 877 879 CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation 878 880 ENDIF 879 CALL histwrite( nid_T, "sossheig", it, ssh n, ndim_hT, ndex_hT ) ! sea surface height881 CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height 880 882 CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux 881 883 CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs … … 884 886 ! in linear free surface case) 885 887 IF( ln_linssh ) THEN 886 zw2d(:,:) = emp (:,:) * ts n(:,:,1,jp_tem)888 zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm) 887 889 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst 888 zw2d(:,:) = emp (:,:) * ts n(:,:,1,jp_sal)890 zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm) 889 891 CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss 890 892 ENDIF … … 922 924 CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping 923 925 CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping 924 IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts n(:,:,1,jp_sal) * tmask(:,:,1)926 IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 925 927 CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping 926 928 ENDIF … … 928 930 CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping 929 931 CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping 930 IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts n(:,:,1,jp_sal) * tmask(:,:,1)932 IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 931 933 CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping 932 934 ENDIF … … 941 943 #endif 942 944 943 CALL histwrite( nid_U, "vozocrtx", it, u n, ndim_U , ndex_U ) ! i-current945 CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current 944 946 CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress 945 947 946 CALL histwrite( nid_V, "vomecrty", it, v n, ndim_V , ndex_V ) ! j-current948 CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current 947 949 CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress 948 950 949 CALL histwrite( nid_W, "vovecrtz", it, w n, ndim_T, ndex_T ) ! vert. current951 CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current 950 952 CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef. 951 953 CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef. … … 974 976 #endif 975 977 976 SUBROUTINE dia_wri_state( cdfile_name )978 SUBROUTINE dia_wri_state( Kmm, cdfile_name ) 977 979 !!--------------------------------------------------------------------- 978 980 !! *** ROUTINE dia_wri_state *** … … 987 989 !! File 'output.abort.nc' is created in case of abnormal job end 988 990 !!---------------------------------------------------------------------- 991 INTEGER , INTENT( in ) :: Kmm ! time level index 989 992 CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created 990 993 !! … … 1003 1006 #endif 1004 1007 1005 CALL iom_rstput( 0, 0, inum, 'votemper', ts n(:,:,:,jp_tem) ) ! now temperature1006 CALL iom_rstput( 0, 0, inum, 'vosaline', ts n(:,:,:,jp_sal) ) ! now salinity1007 CALL iom_rstput( 0, 0, inum, 'sossheig', ssh n) ! sea surface height1008 CALL iom_rstput( 0, 0, inum, 'vozocrtx', u n) ! now i-velocity1009 CALL iom_rstput( 0, 0, inum, 'vomecrty', v n) ! now j-velocity1010 CALL iom_rstput( 0, 0, inum, 'vovecrtz', w n) ! now k-velocity1008 CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature 1009 CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity 1010 CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm) ) ! sea surface height 1011 CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity 1012 CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity 1013 CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity 1011 1014 IF( ALLOCATED(ahtu) ) THEN 1012 1015 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point … … 1024 1027 CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress 1025 1028 IF( .NOT.ln_linssh ) THEN 1026 CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept _n) ! T-cell depth1027 CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t _n) ! T-cell thickness1029 CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm) ) ! T-cell depth 1030 CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm) ) ! T-cell thickness 1028 1031 END IF 1029 1032 IF( ln_wave .AND. ln_sdw ) THEN -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90
r10425 r11771 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 10 11 !!---------------------------------------------------------------------- 11 12 … … 13 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 14 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 !! dom_vvl_sf_ swp: Swap vertical scale factors and update the vertical grid16 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 16 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 17 18 !! dom_vvl_rst : read/write restart file … … 37 38 PUBLIC dom_vvl_init ! called by domain.F90 38 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_ swp! called by step.F9040 PUBLIC dom_vvl_sf_update ! called by step.F90 40 41 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 42 … … 93 94 94 95 95 SUBROUTINE dom_vvl_init 96 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 97 !!---------------------------------------------------------------------- 97 98 !! *** ROUTINE dom_vvl_init *** … … 104 105 !! 105 106 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 106 !! - Regrid: e3 (u/v)_n107 !! e3 (u/v)_b108 !! e3w _n109 !! e3 (u/v)w_b110 !! e3 (u/v)w_n111 !! gdept _n, gdepw_n and gde3w_n107 !! - Regrid: e3[u/v](:,:,:,Kmm) 108 !! e3[u/v](:,:,:,Kmm) 109 !! e3w(:,:,:,Kmm) 110 !! e3[u/v]w_b 111 !! e3[u/v]w_n 112 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 112 113 !! - h(t/u/v)_0 113 114 !! - frq_rst_e3t and frq_rst_hdv … … 115 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 117 !!---------------------------------------------------------------------- 118 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 119 ! 117 120 INTEGER :: ji, jj, jk 118 121 INTEGER :: ii0, ii1, ij0, ij1 … … 130 133 ! 131 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 132 CALL dom_vvl_rst( nit000, 'READ' )133 e3t _a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all135 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 136 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 134 137 ! 135 138 ! !== Set of all other vertical scale factors ==! (now and before) 136 139 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t _b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U138 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )139 CALL dom_vvl_interpol( e3t _b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V140 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )141 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F140 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 141 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 142 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 143 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 144 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 145 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W144 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b (:,:,:), 'W' )145 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW146 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )147 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW148 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )146 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 147 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 148 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 149 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 150 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 151 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 152 150 153 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t _a(:,:,:) = e3t_n(:,:,:)152 e3u _a(:,:,:) = e3u_n(:,:,:)153 e3v _a(:,:,:) = e3v_n(:,:,:)154 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 155 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 156 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 157 ! 155 158 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 156 gdept _n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration)157 gdepw _n(:,:,1) = 0.0_wp158 gde3w _n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg159 gdept _b(:,:,1) = 0.5_wp * e3w_b(:,:,1)160 gdepw _b(:,:,1) = 0.0_wp159 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 160 gdepw(:,:,1,Kmm) = 0.0_wp 161 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 162 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 163 gdepw(:,:,1,Kbb) = 0.0_wp 161 164 DO jk = 2, jpk ! vertical sum 162 165 DO jj = 1,jpj … … 165 168 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 169 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept _n as well as gde3w_ndoes not include the thickness of ISF ??170 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 168 171 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 169 gdepw _n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)170 gdept _n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &171 & + (1-zcoef) * ( gdept _n(ji,jj,jk-1) + e3w_n(ji,jj,jk))172 gde3w _n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)173 gdepw _b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)174 gdept _b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &175 & + (1-zcoef) * ( gdept _b(ji,jj,jk-1) + e3w_b(ji,jj,jk))172 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 173 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 174 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 175 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 176 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 177 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 178 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 176 179 END DO 177 180 END DO … … 179 182 ! 180 183 ! !== thickness of the water column !! (ocean portion only) 181 ht _n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....182 hu _b(:,:) = e3u_b(:,:,1) * umask(:,:,1)183 hu _n(:,:) = e3u_n(:,:,1) * umask(:,:,1)184 hv _b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)185 hv _n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)184 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 185 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 186 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 187 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 188 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 186 189 DO jk = 2, jpkm1 187 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)188 hu _b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)189 hu _n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)190 hv _b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)191 hv _n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)190 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 191 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 192 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 193 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 194 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 192 195 END DO 193 196 ! 194 197 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF196 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )197 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )198 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )198 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 199 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 200 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 201 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 199 202 200 203 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 266 269 267 270 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )271 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 272 !!---------------------------------------------------------------------- 270 273 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 291 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 292 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 293 INTEGER, INTENT( in ) :: kt ! time step 294 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 295 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 296 ! 293 297 INTEGER :: ji, jj, jk ! dummy loop indices … … 321 325 ! ! --------------------------------------------- ! 322 326 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )327 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 328 DO jk = 1, jpkm1 325 ! formally this is the same as e3t _a= e3t_0*(1+ssha/ht_0)326 e3t _a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)329 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 330 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 331 END DO 328 332 ! … … 337 341 zht(:,:) = 0._wp 338 342 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)343 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 344 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 345 END DO 342 346 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 348 352 DO jk = 1, jpkm1 349 353 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )354 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 355 END DO 352 356 ENDIF … … 361 365 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 366 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 368 END DO 365 369 ELSE ! layer case 366 370 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)371 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 372 END DO 369 373 ENDIF … … 476 480 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 481 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )482 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 483 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)484 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 485 END DO 482 486 … … 486 490 ! ! ---baroclinic part--------- ! 487 491 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)492 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 493 END DO 490 494 ENDIF … … 501 505 zht(:,:) = 0.0_wp 502 506 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)504 END DO 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh n(:,:) - zht(:,:) ) )507 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 508 END DO 509 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 510 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t _n))) =', z_tmax511 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 512 ! 509 513 zht(:,:) = 0.0_wp 510 514 DO jk = 1, jpkm1 511 zht(:,:) = zht(:,:) + e3t _a(:,:,jk) * tmask(:,:,jk)512 END DO 513 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh a(:,:) - zht(:,:) ) )515 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 516 END DO 517 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 518 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 515 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t _a))) =', z_tmax519 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 520 ! 517 521 zht(:,:) = 0.0_wp 518 522 DO jk = 1, jpkm1 519 zht(:,:) = zht(:,:) + e3t _b(:,:,jk) * tmask(:,:,jk)520 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh b(:,:) - zht(:,:) ) )523 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 524 END DO 525 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 526 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 523 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t _b))) =', z_tmax524 ! 525 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh b(:,:) ) )527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 530 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh b))) =', z_tmax528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh n(:,:) ) )531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 534 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh n))) =', z_tmax532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh a(:,:) ) )535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 536 ! 537 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 538 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax539 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 540 END IF 537 541 … … 540 544 ! *********************************** ! 541 545 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )546 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 547 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 548 545 549 ! *********************************** ! … … 547 551 ! *********************************** ! 548 552 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)553 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 554 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 555 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)556 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 557 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 558 END DO 555 559 ! ! Inverse of the local depth 556 560 !!gm BUG ? don't understand the use of umask_i here ..... 557 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )558 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )561 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 562 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 563 ! 560 564 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 567 564 568 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***569 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 570 !!---------------------------------------------------------------------- 571 !! *** ROUTINE dom_vvl_sf_update *** 568 572 !! 569 !! ** Purpose : compute time filter and swap of scale factors573 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 574 !! compute all depths and related variables for next time step 571 575 !! write outputs and restart file 572 576 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 577 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 578 !! - reconstruct scale factor at other grid points (interpolate) 575 579 !! - recompute depths and water height fields 576 580 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step581 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 582 !! - Recompute: 579 583 !! e3(u/v)_b 580 !! e3w _n584 !! e3w(:,:,:,Kmm) 581 585 !! e3(u/v)w_b 582 586 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n587 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 588 !! h(u/v) and h(u/v)r 585 589 !! … … 587 591 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 592 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 593 INTEGER, INTENT( in ) :: kt ! time step 594 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 595 ! 591 596 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 600 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 601 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')602 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 603 ! 599 604 IF( kt == nit000 ) THEN 600 605 IF(lwp) WRITE(numout,*) 601 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'602 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'606 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 607 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 608 ENDIF 604 609 ! … … 615 620 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 621 ENDIF 617 gdept _b(:,:,:) = gdept_n(:,:,:)618 gdepw _b(:,:,:) = gdepw_n(:,:,:)619 620 e3t _n(:,:,:) = e3t_a(:,:,:)621 e3u _n(:,:,:) = e3u_a(:,:,:)622 e3v _n(:,:,:) = e3v_a(:,:,:)622 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 623 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 624 625 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 626 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 627 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 623 628 624 629 ! Compute all missing vertical scale factor and depths … … 626 631 ! Horizontal scale factor interpolations 627 632 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also633 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 634 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 635 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )636 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 637 633 638 ! Vertical scale factor interpolations 634 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n(:,:,:), 'W' )635 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' )636 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' )637 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b(:,:,:), 'W' )638 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )639 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )639 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 640 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 641 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 642 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 643 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 644 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 645 641 646 ! t- and w- points depth (set the isf depth as it is in the initial step) 642 gdept _n(:,:,1) = 0.5_wp * e3w_n(:,:,1)643 gdepw _n(:,:,1) = 0.0_wp644 gde3w _n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)647 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 648 gdepw(:,:,1,Kmm) = 0.0_wp 649 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 645 650 DO jk = 2, jpk 646 651 DO jj = 1,jpj … … 649 654 ! 1 for jk = mikt 650 655 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 651 gdepw _n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)652 gdept _n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) &653 & + (1-zcoef) * ( gdept _n(ji,jj,jk-1) + e3w_n(ji,jj,jk) )654 gde3w _n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)656 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 657 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 658 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 659 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 655 660 END DO 656 661 END DO … … 659 664 ! Local depth and Inverse of the local depth of the water 660 665 ! ------------------------------------------------------- 661 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 662 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 ! 664 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 666 ! 667 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 668 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)669 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 670 END DO 668 671 669 672 ! write restart file 670 673 ! ================== 671 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )672 ! 673 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')674 ! 675 END SUBROUTINE dom_vvl_sf_ swp674 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 675 ! 676 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 677 ! 678 END SUBROUTINE dom_vvl_sf_update 676 679 677 680 … … 783 786 784 787 785 SUBROUTINE dom_vvl_rst( kt, cdrw )788 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 789 !!--------------------------------------------------------------------- 787 790 !! *** ROUTINE dom_vvl_rst *** … … 795 798 !! they are set to 0. 796 799 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 800 INTEGER , INTENT(in) :: kt ! ocean time-step 801 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 802 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 803 ! 800 804 INTEGER :: ji, jj, jk … … 806 810 IF( ln_rstart ) THEN !* Read the restart file 807 811 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh n, ldxios = lrxios )812 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 813 ! 810 814 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 817 821 ! ! --------- ! 818 822 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 819 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )820 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 824 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 825 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'826 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 827 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)828 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 829 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 830 END WHERE 827 831 IF( neuler == 0 ) THEN 828 e3t _b(:,:,:) = e3t_n(:,:,:)832 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 833 ENDIF 830 834 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'835 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 836 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 837 IF(lwp) write(numout,*) 'neuler is forced to 0' 834 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )835 e3t _n(:,:,:) = e3t_b(:,:,:)838 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 839 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 836 840 neuler = 0 837 841 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'842 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 843 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 844 IF(lwp) write(numout,*) 'neuler is forced to 0' 841 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )842 e3t _b(:,:,:) = e3t_n(:,:,:)845 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 846 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 843 847 neuler = 0 844 848 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'849 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 850 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 851 IF(lwp) write(numout,*) 'neuler is forced to 0' 848 852 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &853 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 854 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 855 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 856 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)857 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 854 858 neuler = 0 855 859 ENDIF … … 888 892 IF( cn_cfg == 'wad' ) THEN 889 893 ! Wetting and drying test case 890 CALL usr_def_istate( gdept _b, tmask, tsb, ub, vb, sshb)891 ts n (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones892 ssh n (:,:) = sshb(:,:)893 u n (:,:,:) = ub (:,:,:)894 v n (:,:,:) = vb (:,:,:)894 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 895 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 896 ssh (:,:,Kmm) = ssh(:,:,Kbb) 897 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 898 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 899 ELSE 896 900 ! if not test case 897 ssh n(:,:) = -ssh_ref898 ssh b(:,:) = -ssh_ref901 ssh(:,:,Kmm) = -ssh_ref 902 ssh(:,:,Kbb) = -ssh_ref 899 903 900 904 DO jj = 1, jpj 901 905 DO ji = 1, jpi 902 906 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 903 904 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 906 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 907 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 908 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 907 909 ENDIF 908 910 ENDDO … … 912 914 ! Adjust vertical metrics for all wad 913 915 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &916 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 917 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 918 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 919 END DO 918 e3t _b(:,:,:) = e3t_n(:,:,:)920 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 919 921 920 922 DO ji = 1, jpi … … 928 930 ELSE 929 931 ! 930 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t _b and e3t_n931 CALL usr_def_istate( gdept_0, tmask, ts b, ub, vb, sshb )932 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 933 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 932 934 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 933 935 ! 934 936 DO jk=1,jpk 935 e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &936 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &937 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b != 0 on land points937 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 938 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 939 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b != 0 on land points 938 940 END DO 939 e3t_n(:,:,:) = e3t_b(:,:,:) 940 sshn(:,:) = sshb(:,:) ! needed later for gde3w 941 !!$ e3t_n(:,:,:)=e3t_0(:,:,:) 942 !!$ e3t_b(:,:,:)=e3t_0(:,:,:) 941 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 942 ssh(:,: ,Kmm) = ssh(:,: ,Kbb) ! needed later for gde3w 943 943 ! 944 944 END IF ! end of ll_wd edits … … 958 958 ! ! all cases ! 959 959 ! ! --------- ! 960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t _b(:,:,:), ldxios = lwxios )961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t _n(:,:,:), ldxios = lwxios )960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 962 962 ! ! ----------------------- ! 963 963 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/stpctl.F90
r10572 r11771 42 42 CONTAINS 43 43 44 SUBROUTINE stp_ctl( kt, kindic )44 SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 45 45 !!---------------------------------------------------------------------- 46 46 !! *** ROUTINE stp_ctl *** … … 60 60 !!---------------------------------------------------------------------- 61 61 INTEGER, INTENT(in ) :: kt ! ocean time-step index 62 INTEGER, INTENT(in ) :: Kbb, Kmm ! ocean time level index 62 63 INTEGER, INTENT(inout) :: kindic ! error indicator 63 64 !! … … 111 112 ! !== test of extrema ==! 112 113 IF( ll_wd ) THEN 113 zmax(1) = MAXVAL( ABS( ssh n(:,:) + ssh_ref*tmask(:,:,1) ) ) ! ssh max114 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) ) ) ! ssh max 114 115 ELSE 115 zmax(1) = MAXVAL( ABS( ssh n(:,:) ) ) ! ssh max116 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) ) ) ! ssh max 116 117 ENDIF 117 zmax(2) = MAXVAL( ABS( u n(:,:,:) ) ) ! velocity max (zonal only)118 zmax(3) = MAXVAL( -ts n(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp ) ! minus salinity max119 zmax(4) = MAXVAL( ts n(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp ) ! salinity max120 zmax(5) = MAXVAL( -ts n(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp ) ! minus temperature max121 zmax(6) = MAXVAL( ts n(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp ) ! temperature max118 zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) ) ! velocity max (zonal only) 119 zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus salinity max 120 zmax(4) = MAXVAL( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! salinity max 121 zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus temperature max 122 zmax(6) = MAXVAL( ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! temperature max 122 123 zmax(7) = REAL( nstop , wp ) ! stop indicator 123 124 IF( ln_zad_Aimp ) THEN … … 148 149 ! !== error handling ==! 149 150 IF( ( ln_ctl .OR. lsomeoce ) .AND. ( & ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 150 & zmax(1) > 50._wp .OR. & ! too large sea surface height ( > 50 m )151 & zmax(2) > 20._wp .OR. & ! too large velocity ( > 20 m/s)152 !!$ 153 !!$ 154 !!$ 151 & zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 152 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 153 !!$ & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 154 !!$ & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 155 !!$ & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 155 156 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN ! NaN encounter in the tests 156 157 IF( lk_mpp .AND. ln_ctl ) THEN 157 CALL mpp_maxloc( 'stpctl', ABS(ssh n) , ssmask(:,:) , zzz, ih )158 CALL mpp_maxloc( 'stpctl', ABS(u n) , umask (:,:,:), zzz, iu )159 CALL mpp_minloc( 'stpctl', ts n(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 )160 CALL mpp_maxloc( 'stpctl', ts n(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 )158 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm)) , ssmask(:,:) , zzz, ih ) 159 CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm)) , umask (:,:,:), zzz, iu ) 160 CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 161 CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 161 162 ELSE 162 ih(:) = MAXLOC( ABS( ssh n(:,:) ) ) + (/ nimpp - 1, njmpp - 1 /)163 iu(:) = MAXLOC( ABS( u n (:,:,:) ) ) + (/ nimpp - 1, njmpp - 1, 0 /)164 is1(:) = MINLOC( ts n(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)165 is2(:) = MAXLOC( ts n(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)163 ih(:) = MAXLOC( ABS( ssh(:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1 /) 164 iu(:) = MAXLOC( ABS( uu (:,:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1, 0 /) 165 is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 166 is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 166 167 ENDIF 167 168 168 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 50 m or |U| > 20 m/sor NaN encounter in the tests'169 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 169 170 WRITE(ctmp2,9100) kt, zmax(1), ih(1) , ih(2) 170 171 WRITE(ctmp3,9200) kt, zmax(2), iu(1) , iu(2) , iu(3) … … 173 174 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort.nc file' 174 175 175 CALL dia_wri_state( 'output.abort' ) ! create an output.abort file176 CALL dia_wri_state( Kmm, 'output.abort' ) ! create an output.abort file 176 177 177 178 IF( .NOT. ln_ctl ) THEN -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/trazdf.F90
r10572 r11771 44 44 CONTAINS 45 45 46 SUBROUTINE tra_zdf( kt )46 SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_zdf *** … … 50 50 !! ** Purpose : compute the vertical ocean tracer physics. 51 51 !!--------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER , INTENT(in) :: kt ! ocean time-step index 53 INTEGER , INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices 54 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 53 55 ! 54 56 INTEGER :: jk ! Dummy loop indices … … 70 72 IF( l_trdtra ) THEN !* Save ta and sa trends 71 73 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)74 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 75 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 74 76 ENDIF 75 77 ! 76 78 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )79 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 78 80 79 81 !!gm WHY here ! and I don't like that ! … … 81 83 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 84 ! JMM : restore negative salinities to small salinities: 83 !!$ WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp85 !!$ WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 84 86 !!gm 85 87 86 88 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 89 DO jk = 1, jpkm1 88 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &89 & / (e3t _n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)90 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &91 & / (e3t _n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)90 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 91 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 92 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 92 94 END DO 93 95 !!gm this should be moved in trdtra.F90 and done on all trends 94 96 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. ) 95 97 !!gm 96 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )97 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )98 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 99 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 98 100 DEALLOCATE( ztrdt , ztrds ) 99 101 ENDIF 100 102 ! ! print mean trends (used for debugging) 101 IF(ln_ctl) CALL prt_ctl( tab3d_1= tsa(:,:,:,jp_tem), clinfo1=' zdf - Ta: ', mask1=tmask, &102 & tab3d_2= tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )103 IF(ln_ctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, & 104 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 103 105 ! 104 106 IF( ln_timing ) CALL timing_stop('tra_zdf') … … 107 109 108 110 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )111 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 110 112 !!---------------------------------------------------------------------- 111 113 !! *** ROUTINE tra_zdf_imp *** … … 125 127 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 128 !! 127 !! ** Action : - pt abecomes the after tracer128 !!--------------------------------------------------------------------- 129 INTEGER , INTENT(in ) :: kt ! ocean time-step index130 INTEGER , INTENT(in ) :: kit000 ! first time step index131 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)132 INTEGER , INTENT(in ) :: kjpt ! number of tracers133 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step134 REAL(wp) , DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields135 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt ), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field129 !! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer 130 !!--------------------------------------------------------------------- 131 INTEGER , INTENT(in ) :: kt ! ocean time-step index 132 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 133 INTEGER , INTENT(in ) :: kit000 ! first time step index 134 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 135 INTEGER , INTENT(in ) :: kjpt ! number of tracers 136 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 137 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 136 138 ! 137 139 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 181 183 DO jj = 2, jpjm1 182 184 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 183 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w _n(ji,jj,jk)184 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w _n(ji,jj,jk+1)185 zwd(ji,jj,jk) = e3t _a(ji,jj,jk) - zzwi - zzws &185 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 186 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 187 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 186 188 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 187 189 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) … … 194 196 DO jj = 2, jpjm1 195 197 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w _n(ji,jj,jk)197 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w _n(ji,jj,jk+1)198 zwd(ji,jj,jk) = e3t _a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)198 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 199 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 200 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 199 201 END DO 200 202 END DO … … 218 220 ! The solution will be in the 4d array pta. 219 221 ! The 3d array zwt is used as a work space array. 220 ! En route to the solution pt ais used a to evaluate the rhs and then222 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 221 223 ! used as a work space array: its value is modified. 222 224 ! … … 238 240 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 239 241 DO ji = fs_2, fs_jpim1 240 pt a(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)242 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 241 243 END DO 242 244 END DO … … 244 246 DO jj = 2, jpjm1 245 247 DO ji = fs_2, fs_jpim1 246 zrhs = e3t _b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn) ! zrhs=right hand side247 pt a(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)248 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 249 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 248 250 END DO 249 251 END DO … … 252 254 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 253 255 DO ji = fs_2, fs_jpim1 254 pt a(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)256 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 255 257 END DO 256 258 END DO … … 258 260 DO jj = 2, jpjm1 259 261 DO ji = fs_2, fs_jpim1 260 pt a(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &262 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 261 263 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 262 264 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/usrdef_sbc.F90
r10976 r11771 40 40 CONTAINS 41 41 42 SUBROUTINE usrdef_sbc_oce( kt, K bb )42 SUBROUTINE usrdef_sbc_oce( kt, Kmm, Kbb ) 43 43 !!--------------------------------------------------------------------- 44 44 !! *** ROUTINE usr_def_sbc *** … … 54 54 !! 55 55 !!---------------------------------------------------------------------- 56 INTEGER, INTENT(in) :: kt ! ocean time step57 INTEGER, INTENT(in) :: Kbb ! ocean time index56 INTEGER, INTENT(in) :: kt ! ocean time step 57 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time index 58 58 INTEGER :: ji, jj ! dummy loop indices 59 59 REAL(wp) :: zrhoair = 1.22 ! approximate air density [Kg/m3] … … 88 88 89 89 WHERE( ABS(gphit) <= rn_windszy/2. ) 90 zwndrel(:,:) = rn_u10 - rn_uofac * u n(:,:,1)90 zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kmm) 91 91 ELSEWHERE 92 zwndrel(:,:) = - rn_uofac * u n(:,:,1)92 zwndrel(:,:) = - rn_uofac * uu(:,:,1,Kmm) 93 93 END WHERE 94 94 utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 95 95 96 zwndrel(:,:) = - rn_uofac * v n(:,:,1)96 zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kmm) 97 97 vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 98 98
Note: See TracChangeset
for help on using the changeset viewer.
