Changeset 12377 for NEMO/trunk/src/OCE/DYN/dynhpg.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynhpg.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/DYN/dynhpg.F90
r11536 r12377 31 31 !!---------------------------------------------------------------------- 32 32 USE oce ! ocean dynamics and tracers 33 USE isf_oce , ONLY : risfload ! ice shelf (risfload variable) 34 USE isfload , ONLY : isf_load ! ice shelf (isf_load routine ) 33 35 USE sbc_oce ! surface variable (only for the flag with ice shelf) 34 36 USE dom_oce ! ocean space and time domain … … 73 75 74 76 !! * Substitutions 75 # include " vectopt_loop_substitute.h90"77 # include "do_loop_substitute.h90" 76 78 !!---------------------------------------------------------------------- 77 79 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 81 83 CONTAINS 82 84 83 SUBROUTINE dyn_hpg( kt )85 SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) 84 86 !!--------------------------------------------------------------------- 85 87 !! *** ROUTINE dyn_hpg *** … … 88 90 !! using the scheme defined in the namelist 89 91 !! 90 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend92 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 91 93 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 92 94 !!---------------------------------------------------------------------- 93 INTEGER, INTENT(in) :: kt ! ocean time-step index 95 INTEGER , INTENT( in ) :: kt ! ocean time-step index 96 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 97 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 98 ! 94 99 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 95 100 !!---------------------------------------------------------------------- … … 97 102 IF( ln_timing ) CALL timing_start('dyn_hpg') 98 103 ! 99 IF( l_trddyn ) THEN ! Temporary saving of ua and vatrends (l_trddyn)104 IF( l_trddyn ) THEN ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) 100 105 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 101 ztrdu(:,:,:) = ua(:,:,:)102 ztrdv(:,:,:) = va(:,:,:)106 ztrdu(:,:,:) = puu(:,:,:,Krhs) 107 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 103 108 ENDIF 104 109 ! 105 110 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 106 CASE ( np_zco ) ; CALL hpg_zco ( kt )! z-coordinate107 CASE ( np_zps ) ; CALL hpg_zps ( kt )! z-coordinate plus partial steps (interpolation)108 CASE ( np_sco ) ; CALL hpg_sco ( kt )! s-coordinate (standard jacobian formulation)109 CASE ( np_djc ) ; CALL hpg_djc ( kt )! s-coordinate (Density Jacobian with Cubic polynomial)110 CASE ( np_prj ) ; CALL hpg_prj ( kt )! s-coordinate (Pressure Jacobian scheme)111 CASE ( np_isf ) ; CALL hpg_isf ( kt )! s-coordinate similar to sco modify for ice shelf111 CASE ( np_zco ) ; CALL hpg_zco ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate 112 CASE ( np_zps ) ; CALL hpg_zps ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate plus partial steps (interpolation) 113 CASE ( np_sco ) ; CALL hpg_sco ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (standard jacobian formulation) 114 CASE ( np_djc ) ; CALL hpg_djc ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial) 115 CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme) 116 CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf 112 117 END SELECT 113 118 ! 114 119 IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 115 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)116 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)117 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )120 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 121 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 122 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt, Kmm ) 118 123 DEALLOCATE( ztrdu , ztrdv ) 119 124 ENDIF 120 125 ! 121 IF( ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, &122 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )126 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg - Ua: ', mask1=umask, & 127 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 123 128 ! 124 129 IF( ln_timing ) CALL timing_stop('dyn_hpg') … … 127 132 128 133 129 SUBROUTINE dyn_hpg_init 134 SUBROUTINE dyn_hpg_init( Kmm ) 130 135 !!---------------------------------------------------------------------- 131 136 !! *** ROUTINE dyn_hpg_init *** … … 137 142 !! with the type of vertical coordinate used (zco, zps, sco) 138 143 !!---------------------------------------------------------------------- 144 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 145 ! 139 146 INTEGER :: ioptio = 0 ! temporary integer 140 147 INTEGER :: ios ! Local integer output status for namelist read … … 150 157 !!---------------------------------------------------------------------- 151 158 ! 152 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient153 159 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 154 160 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' ) 155 161 ! 156 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient157 162 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 158 163 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' ) … … 213 218 ENDIF 214 219 ! 215 IF ( .NOT. ln_isfcav ) THEN !--- no ice shelf load216 riceload(:,:) = 0._wp217 !218 ELSE !--- set an ice shelf load219 !220 IF(lwp) WRITE(numout,*)221 IF(lwp) WRITE(numout,*) ' ice shelf case: set the ice-shelf load'222 ALLOCATE( zts_top(jpi,jpj,jpts) , zrhd(jpi,jpj,jpk) , zrhdtop_isf(jpi,jpj) , ziceload(jpi,jpj) )223 !224 znad = 1._wp !- To use density and not density anomaly225 !226 ! !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)227 zts_top(:,:,jp_tem) = -1.9_wp ; zts_top(:,:,jp_sal) = 34.4_wp228 !229 DO jk = 1, jpk !- compute density of the water displaced by the ice shelf230 CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )231 END DO232 !233 ! !- compute rhd at the ice/oce interface (ice shelf side)234 CALL eos( zts_top , risfdep, zrhdtop_isf )235 !236 ! !- Surface value + ice shelf gradient237 ziceload = 0._wp ! compute pressure due to ice shelf load238 DO jj = 1, jpj ! (used to compute hpgi/j for all the level from 1 to miku/v)239 DO ji = 1, jpi ! divided by 2 later240 ikt = mikt(ji,jj)241 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))242 DO jk = 2, ikt-1243 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &244 & * (1._wp - tmask(ji,jj,jk))245 END DO246 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &247 & * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )248 END DO249 END DO250 riceload(:,:) = ziceload(:,:) ! need to be saved for diaar5251 !252 DEALLOCATE( zts_top , zrhd , zrhdtop_isf , ziceload )253 ENDIF254 !255 220 END SUBROUTINE dyn_hpg_init 256 221 257 222 258 SUBROUTINE hpg_zco( kt )223 SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) 259 224 !!--------------------------------------------------------------------- 260 225 !! *** ROUTINE hpg_zco *** … … 266 231 !! level: zhpi = grav ..... 267 232 !! zhpj = grav ..... 268 !! add it to the general momentum trend (ua,va). 269 !! ua = ua - 1/e1u * zhpi 270 !! va = va - 1/e2v * zhpj 271 !! 272 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 273 !!---------------------------------------------------------------------- 274 INTEGER, INTENT(in) :: kt ! ocean time-step index 233 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 234 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 235 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 236 !! 237 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 238 !!---------------------------------------------------------------------- 239 INTEGER , INTENT( in ) :: kt ! ocean time-step index 240 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 241 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 275 242 ! 276 243 INTEGER :: ji, jj, jk ! dummy loop indices … … 288 255 289 256 ! Surface value 290 DO jj = 2, jpjm1 291 DO ji = fs_2, fs_jpim1 ! vector opt. 292 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 293 ! hydrostatic pressure gradient 294 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 295 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 296 ! add to the general momentum trend 297 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 298 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 299 END DO 300 END DO 257 DO_2D_00_00 258 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 259 ! hydrostatic pressure gradient 260 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 261 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 262 ! add to the general momentum trend 263 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 264 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 265 END_2D 301 266 302 267 ! 303 268 ! interior value (2=<jk=<jpkm1) 304 DO jk = 2, jpkm1 305 DO jj = 2, jpjm1 306 DO ji = fs_2, fs_jpim1 ! vector opt. 307 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 308 ! hydrostatic pressure gradient 309 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 310 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 311 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 312 313 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 314 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 315 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 316 ! add to the general momentum trend 317 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 318 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 319 END DO 320 END DO 321 END DO 269 DO_3D_00_00( 2, jpkm1 ) 270 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 271 ! hydrostatic pressure gradient 272 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 273 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 274 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 275 276 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 277 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 278 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 279 ! add to the general momentum trend 280 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 281 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 282 END_3D 322 283 ! 323 284 END SUBROUTINE hpg_zco 324 285 325 286 326 SUBROUTINE hpg_zps( kt )287 SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs ) 327 288 !!--------------------------------------------------------------------- 328 289 !! *** ROUTINE hpg_zps *** … … 330 291 !! ** Method : z-coordinate plus partial steps case. blahblah... 331 292 !! 332 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 333 !!---------------------------------------------------------------------- 334 INTEGER, INTENT(in) :: kt ! ocean time-step index 293 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 294 !!---------------------------------------------------------------------- 295 INTEGER , INTENT( in ) :: kt ! ocean time-step index 296 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 297 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 335 298 !! 336 299 INTEGER :: ji, jj, jk ! dummy loop indices … … 348 311 349 312 ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level 350 CALL zps_hde( kt, jpts, tsn, zgtsu, zgtsv, rhd, zgru , zgrv )313 CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv ) 351 314 352 315 ! Local constant initialization … … 354 317 355 318 ! Surface value (also valid in partial step case) 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 ! vector opt. 358 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 359 ! hydrostatic pressure gradient 360 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 361 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 362 ! add to the general momentum trend 363 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 364 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 365 END DO 366 END DO 319 DO_2D_00_00 320 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 321 ! hydrostatic pressure gradient 322 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 323 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 324 ! add to the general momentum trend 325 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 326 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 327 END_2D 367 328 368 329 ! interior value (2=<jk=<jpkm1) 369 DO jk = 2, jpkm1 370 DO jj = 2, jpjm1 371 DO ji = fs_2, fs_jpim1 ! vector opt. 372 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 373 ! hydrostatic pressure gradient 374 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 375 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 376 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 377 378 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 379 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 380 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 381 ! add to the general momentum trend 382 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 383 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 384 END DO 385 END DO 386 END DO 330 DO_3D_00_00( 2, jpkm1 ) 331 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 332 ! hydrostatic pressure gradient 333 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 334 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 335 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 336 337 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 338 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 339 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 340 ! add to the general momentum trend 341 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 342 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 343 END_3D 387 344 388 345 ! partial steps correction at the last level (use zgru & zgrv computed in zpshde.F90) 389 DO jj = 2, jpjm1 390 DO ji = 2, jpim1 391 iku = mbku(ji,jj) 392 ikv = mbkv(ji,jj) 393 zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj ,iku) ) 394 zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) ) 395 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 396 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 397 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 398 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 399 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 400 ENDIF 401 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 402 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 403 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 404 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 405 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 406 ENDIF 407 END DO 408 END DO 346 DO_2D_00_00 347 iku = mbku(ji,jj) 348 ikv = mbkv(ji,jj) 349 zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj ,iku,Kmm) ) 350 zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji ,jj+1,ikv,Kmm) ) 351 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 352 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value 353 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 354 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 355 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 356 ENDIF 357 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 358 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value 359 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 360 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 361 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 362 ENDIF 363 END_2D 409 364 ! 410 365 END SUBROUTINE hpg_zps 411 366 412 367 413 SUBROUTINE hpg_sco( kt )368 SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 414 369 !!--------------------------------------------------------------------- 415 370 !! *** ROUTINE hpg_sco *** … … 423 378 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 424 379 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 425 !! add it to the general momentum trend (ua,va). 426 !! ua = ua - 1/e1u * zhpi 427 !! va = va - 1/e2v * zhpj 428 !! 429 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 430 !!---------------------------------------------------------------------- 431 INTEGER, INTENT(in) :: kt ! ocean time-step index 380 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 381 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 382 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 383 !! 384 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 385 !!---------------------------------------------------------------------- 386 INTEGER , INTENT( in ) :: kt ! ocean time-step index 387 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 388 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 432 389 !! 433 390 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices … … 452 409 ! 453 410 IF( ln_wd_il ) THEN 454 DO jj = 2, jpjm1 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 461 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 462 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 464 IF(ll_tmp1) THEN 465 zcpx(ji,jj) = 1.0_wp 466 ELSE IF(ll_tmp2) THEN 467 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 469 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 470 ELSE 471 zcpx(ji,jj) = 0._wp 472 END IF 473 474 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 476 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 477 & > rn_wdmin1 + rn_wdmin2 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 479 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 480 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 482 IF(ll_tmp1) THEN 483 zcpy(ji,jj) = 1.0_wp 484 ELSE IF(ll_tmp2) THEN 485 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 487 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 488 ELSE 489 zcpy(ji,jj) = 0._wp 490 END IF 491 END DO 492 END DO 411 DO_2D_00_00 412 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 413 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 414 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 415 & > rn_wdmin1 + rn_wdmin2 416 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 417 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 418 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 419 420 IF(ll_tmp1) THEN 421 zcpx(ji,jj) = 1.0_wp 422 ELSE IF(ll_tmp2) THEN 423 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 424 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 425 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 426 ELSE 427 zcpx(ji,jj) = 0._wp 428 END IF 429 430 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 431 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 432 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 433 & > rn_wdmin1 + rn_wdmin2 434 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 435 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 436 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 437 438 IF(ll_tmp1) THEN 439 zcpy(ji,jj) = 1.0_wp 440 ELSE IF(ll_tmp2) THEN 441 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 442 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 443 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 444 ELSE 445 zcpy(ji,jj) = 0._wp 446 END IF 447 END_2D 493 448 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. ) 494 449 END IF 495 450 496 451 ! Surface value 497 DO jj = 2, jpjm1 498 DO ji = fs_2, fs_jpim1 ! vector opt. 499 ! hydrostatic pressure gradient along s-surfaces 500 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 501 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 502 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 503 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 504 ! s-coordinate pressure gradient correction 505 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 506 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 507 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 508 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 509 ! 510 IF( ln_wd_il ) THEN 511 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 512 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 513 zuap = zuap * zcpx(ji,jj) 514 zvap = zvap * zcpy(ji,jj) 515 ENDIF 516 ! 517 ! add to the general momentum trend 518 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 519 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 520 END DO 521 END DO 452 DO_2D_00_00 453 ! hydrostatic pressure gradient along s-surfaces 454 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 455 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 456 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 457 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 458 ! s-coordinate pressure gradient correction 459 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 460 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 461 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 462 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 463 ! 464 IF( ln_wd_il ) THEN 465 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 466 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 467 zuap = zuap * zcpx(ji,jj) 468 zvap = zvap * zcpy(ji,jj) 469 ENDIF 470 ! 471 ! add to the general momentum trend 472 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 473 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 474 END_2D 522 475 523 476 ! interior value (2=<jk=<jpkm1) 524 DO jk = 2, jpkm1 525 DO jj = 2, jpjm1 526 DO ji = fs_2, fs_jpim1 ! vector opt. 527 ! hydrostatic pressure gradient along s-surfaces 528 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 529 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 530 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 531 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 532 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 533 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 534 ! s-coordinate pressure gradient correction 535 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 536 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 537 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 538 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 539 ! 540 IF( ln_wd_il ) THEN 541 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 542 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 543 zuap = zuap * zcpx(ji,jj) 544 zvap = zvap * zcpy(ji,jj) 545 ENDIF 546 ! 547 ! add to the general momentum trend 548 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 549 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 550 END DO 551 END DO 552 END DO 477 DO_3D_00_00( 2, jpkm1 ) 478 ! hydrostatic pressure gradient along s-surfaces 479 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 480 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 481 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 482 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 483 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 484 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 485 ! s-coordinate pressure gradient correction 486 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 487 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 488 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 489 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 490 ! 491 IF( ln_wd_il ) THEN 492 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 493 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 494 zuap = zuap * zcpx(ji,jj) 495 zvap = zvap * zcpy(ji,jj) 496 ENDIF 497 ! 498 ! add to the general momentum trend 499 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 500 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 501 END_3D 553 502 ! 554 503 IF( ln_wd_il ) DEALLOCATE( zcpx , zcpy ) … … 557 506 558 507 559 SUBROUTINE hpg_isf( kt )508 SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 560 509 !!--------------------------------------------------------------------- 561 510 !! *** ROUTINE hpg_isf *** … … 569 518 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 570 519 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 571 !! add it to the general momentum trend ( ua,va).572 !! ua = ua- 1/e1u * zhpi573 !! va = va- 1/e2v * zhpj574 !! iceload is added and partial cell case are added to the top and bottom520 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 521 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 522 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 523 !! iceload is added 575 524 !! 576 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 577 !!---------------------------------------------------------------------- 578 INTEGER, INTENT(in) :: kt ! ocean time-step index 525 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 526 !!---------------------------------------------------------------------- 527 INTEGER , INTENT( in ) :: kt ! ocean time-step index 528 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 529 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 579 530 !! 580 531 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices … … 597 548 DO jj = 1, jpj 598 549 ikt = mikt(ji,jj) 599 zts_top(ji,jj,1) = ts n(ji,jj,ikt,1)600 zts_top(ji,jj,2) = ts n(ji,jj,ikt,2)550 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 551 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 601 552 END DO 602 553 END DO … … 606 557 !===== Compute surface value ===================================================== 607 558 !================================================================================== 608 DO jj = 2, jpjm1 609 DO ji = fs_2, fs_jpim1 ! vector opt. 610 ikt = mikt(ji,jj) 611 iktp1i = mikt(ji+1,jj) 612 iktp1j = mikt(ji,jj+1) 613 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 614 ! we assume ISF is in isostatic equilibrium 615 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i) & 616 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 617 & - 0.5_wp * e3w_n(ji,jj,ikt) & 618 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 619 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 620 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j) & 621 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 622 & - 0.5_wp * e3w_n(ji,jj,ikt) & 623 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 624 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 625 ! s-coordinate pressure gradient correction (=0 if z coordinate) 626 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 627 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 628 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 629 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 630 ! add to the general momentum trend 631 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 632 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 633 END DO 634 END DO 559 DO_2D_00_00 560 ikt = mikt(ji,jj) 561 iktp1i = mikt(ji+1,jj) 562 iktp1j = mikt(ji,jj+1) 563 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 564 ! we assume ISF is in isostatic equilibrium 565 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 566 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 567 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 568 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 569 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 570 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 571 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 572 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 573 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 574 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 575 ! s-coordinate pressure gradient correction (=0 if z coordinate) 576 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 577 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 578 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 579 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 580 ! add to the general momentum trend 581 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 582 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 583 END_2D 635 584 !================================================================================== 636 585 !===== Compute interior value ===================================================== 637 586 !================================================================================== 638 587 ! interior value (2=<jk=<jpkm1) 639 DO jk = 2, jpkm1 640 DO jj = 2, jpjm1 641 DO ji = fs_2, fs_jpim1 ! vector opt. 642 ! hydrostatic pressure gradient along s-surfaces 643 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 644 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 645 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 646 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 647 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 648 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 649 ! s-coordinate pressure gradient correction 650 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 651 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 652 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 653 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 654 ! add to the general momentum trend 655 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 656 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 657 END DO 658 END DO 659 END DO 588 DO_3D_00_00( 2, jpkm1 ) 589 ! hydrostatic pressure gradient along s-surfaces 590 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 591 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 592 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 593 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 594 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 595 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 596 ! s-coordinate pressure gradient correction 597 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 598 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 599 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 600 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 601 ! add to the general momentum trend 602 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 603 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 604 END_3D 660 605 ! 661 606 END SUBROUTINE hpg_isf 662 607 663 608 664 SUBROUTINE hpg_djc( kt )609 SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 665 610 !!--------------------------------------------------------------------- 666 611 !! *** ROUTINE hpg_djc *** … … 670 615 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 671 616 !!---------------------------------------------------------------------- 672 INTEGER, INTENT(in) :: kt ! ocean time-step index 617 INTEGER , INTENT( in ) :: kt ! ocean time-step index 618 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 619 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 673 620 !! 674 621 INTEGER :: ji, jj, jk ! dummy loop indices … … 686 633 IF( ln_wd_il ) THEN 687 634 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 688 DO jj = 2, jpjm1 689 DO ji = 2, jpim1 690 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 691 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 692 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 693 & > rn_wdmin1 + rn_wdmin2 694 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 695 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 696 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 697 IF(ll_tmp1) THEN 698 zcpx(ji,jj) = 1.0_wp 699 ELSE IF(ll_tmp2) THEN 700 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 701 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 702 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 703 ELSE 704 zcpx(ji,jj) = 0._wp 705 END IF 706 707 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 708 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 709 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 710 & > rn_wdmin1 + rn_wdmin2 711 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 712 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 713 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 714 715 IF(ll_tmp1) THEN 716 zcpy(ji,jj) = 1.0_wp 717 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 719 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 720 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 721 ELSE 722 zcpy(ji,jj) = 0._wp 723 END IF 724 END DO 725 END DO 635 DO_2D_00_00 636 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 637 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 638 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 639 & > rn_wdmin1 + rn_wdmin2 640 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 641 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 642 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 643 IF(ll_tmp1) THEN 644 zcpx(ji,jj) = 1.0_wp 645 ELSE IF(ll_tmp2) THEN 646 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 647 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 648 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 649 ELSE 650 zcpx(ji,jj) = 0._wp 651 END IF 652 653 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 654 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 655 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 656 & > rn_wdmin1 + rn_wdmin2 657 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 658 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 659 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 660 661 IF(ll_tmp1) THEN 662 zcpy(ji,jj) = 1.0_wp 663 ELSE IF(ll_tmp2) THEN 664 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 665 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 666 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 667 ELSE 668 zcpy(ji,jj) = 0._wp 669 END IF 670 END_2D 726 671 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. ) 727 672 END IF … … 744 689 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 745 690 746 DO jk = 2, jpkm1 747 DO jj = 2, jpjm1 748 DO ji = fs_2, fs_jpim1 ! vector opt. 749 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 750 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 751 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 752 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 753 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 754 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 755 END DO 756 END DO 757 END DO 691 DO_3D_00_00( 2, jpkm1 ) 692 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 693 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 694 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 695 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 696 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 697 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 698 END_3D 758 699 759 700 !------------------------------------------------------------------------- … … 765 706 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 766 707 767 DO jk = 2, jpkm1 768 DO jj = 2, jpjm1 769 DO ji = fs_2, fs_jpim1 ! vector opt. 770 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 771 772 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 773 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 774 775 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 776 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 777 778 IF( cffw > zep) THEN 779 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 780 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 781 ELSE 782 drhow(ji,jj,jk) = 0._wp 783 ENDIF 784 785 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 786 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 787 788 IF( cffu > zep ) THEN 789 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 790 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 791 ELSE 792 drhou(ji,jj,jk ) = 0._wp 793 ENDIF 794 795 IF( cffx > zep ) THEN 796 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) & 797 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 798 ELSE 799 dzu(ji,jj,jk) = 0._wp 800 ENDIF 801 802 IF( cffv > zep ) THEN 803 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 804 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 805 ELSE 806 drhov(ji,jj,jk) = 0._wp 807 ENDIF 808 809 IF( cffy > zep ) THEN 810 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 811 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 812 ELSE 813 dzv(ji,jj,jk) = 0._wp 814 ENDIF 815 816 END DO 817 END DO 818 END DO 708 DO_3D_00_00( 2, jpkm1 ) 709 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 710 711 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 712 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 713 714 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 715 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 716 717 IF( cffw > zep) THEN 718 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 719 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 720 ELSE 721 drhow(ji,jj,jk) = 0._wp 722 ENDIF 723 724 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 725 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 726 727 IF( cffu > zep ) THEN 728 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 729 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 730 ELSE 731 drhou(ji,jj,jk ) = 0._wp 732 ENDIF 733 734 IF( cffx > zep ) THEN 735 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) & 736 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 737 ELSE 738 dzu(ji,jj,jk) = 0._wp 739 ENDIF 740 741 IF( cffv > zep ) THEN 742 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 743 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 744 ELSE 745 drhov(ji,jj,jk) = 0._wp 746 ENDIF 747 748 IF( cffy > zep ) THEN 749 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 750 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 751 ELSE 752 dzv(ji,jj,jk) = 0._wp 753 ENDIF 754 755 END_3D 819 756 820 757 !---------------------------------------------------------------------------------- … … 837 774 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 838 775 839 DO jj = 2, jpjm1 840 DO ji = fs_2, fs_jpim1 ! vector opt. 841 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) & 842 & * ( rhd(ji,jj,1) & 843 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 844 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) & 845 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) ) 846 END DO 847 END DO 776 DO_2D_00_00 777 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 778 & * ( rhd(ji,jj,1) & 779 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 780 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 781 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 782 END_2D 848 783 849 784 !!bug gm : here also, simplification is possible 850 785 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop 851 786 852 DO jk = 2, jpkm1 853 DO jj = 2, jpjm1 854 DO ji = fs_2, fs_jpim1 ! vector opt. 855 856 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 857 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) & 858 & - grav * z1_10 * ( & 859 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 860 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 861 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 862 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 863 & ) 864 865 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 866 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) & 867 & - grav* z1_10 * ( & 868 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 869 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 870 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 871 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 872 & ) 873 874 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 875 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) & 876 & - grav* z1_10 * ( & 877 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 878 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 879 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 880 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 881 & ) 882 883 END DO 884 END DO 885 END DO 787 DO_3D_00_00( 2, jpkm1 ) 788 789 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 790 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 791 & - grav * z1_10 * ( & 792 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 793 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 794 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 795 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 796 & ) 797 798 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 799 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 800 & - grav* z1_10 * ( & 801 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 802 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 803 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 804 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 805 & ) 806 807 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 808 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 809 & - grav* z1_10 * ( & 810 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 811 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 812 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 813 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 814 & ) 815 816 END_3D 886 817 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1., rho_i, 'U', 1., rho_j, 'V', 1. ) 887 818 … … 889 820 ! Surface value 890 821 ! --------------- 891 DO jj = 2, jpjm1 892 DO ji = fs_2, fs_jpim1 ! vector opt. 893 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 894 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 895 IF( ln_wd_il ) THEN 896 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 897 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 898 ENDIF 899 ! add to the general momentum trend 900 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 901 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 902 END DO 903 END DO 822 DO_2D_00_00 823 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 824 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 825 IF( ln_wd_il ) THEN 826 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 827 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 828 ENDIF 829 ! add to the general momentum trend 830 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 831 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 832 END_2D 904 833 905 834 ! ---------------- 906 835 ! interior value (2=<jk=<jpkm1) 907 836 ! ---------------- 908 DO jk = 2, jpkm1 909 DO jj = 2, jpjm1 910 DO ji = fs_2, fs_jpim1 ! vector opt. 911 ! hydrostatic pressure gradient along s-surfaces 912 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 913 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 914 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 915 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 916 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 917 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 918 IF( ln_wd_il ) THEN 919 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 920 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 921 ENDIF 922 ! add to the general momentum trend 923 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 924 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 925 END DO 926 END DO 927 END DO 837 DO_3D_00_00( 2, jpkm1 ) 838 ! hydrostatic pressure gradient along s-surfaces 839 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 840 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 841 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 842 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 843 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 844 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 845 IF( ln_wd_il ) THEN 846 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 847 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 848 ENDIF 849 ! add to the general momentum trend 850 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 851 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 852 END_3D 928 853 ! 929 854 IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy ) … … 932 857 933 858 934 SUBROUTINE hpg_prj( kt )859 SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 935 860 !!--------------------------------------------------------------------- 936 861 !! *** ROUTINE hpg_prj *** … … 941 866 !! all vertical coordinate systems 942 867 !! 943 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend868 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 944 869 !!---------------------------------------------------------------------- 945 870 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 946 INTEGER, INTENT(in) :: kt ! ocean time-step index 871 INTEGER , INTENT( in ) :: kt ! ocean time-step index 872 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 873 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 947 874 !! 948 875 INTEGER :: ji, jj, jk, jkk ! dummy loop indices … … 974 901 IF( ln_wd_il ) THEN 975 902 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 976 DO jj = 2, jpjm1 977 DO ji = 2, jpim1 978 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 979 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 980 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 981 & > rn_wdmin1 + rn_wdmin2 982 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 983 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 984 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 985 986 IF(ll_tmp1) THEN 987 zcpx(ji,jj) = 1.0_wp 988 ELSE IF(ll_tmp2) THEN 989 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 990 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 991 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 992 993 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 994 ELSE 995 zcpx(ji,jj) = 0._wp 996 END IF 997 998 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 999 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1000 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1001 & > rn_wdmin1 + rn_wdmin2 1002 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 1003 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1004 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1005 1006 IF(ll_tmp1) THEN 1007 zcpy(ji,jj) = 1.0_wp 1008 ELSE IF(ll_tmp2) THEN 1009 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1010 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1011 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1012 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1013 1014 ELSE 1015 zcpy(ji,jj) = 0._wp 1016 ENDIF 1017 END DO 1018 END DO 903 DO_2D_00_00 904 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 905 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 906 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 907 & > rn_wdmin1 + rn_wdmin2 908 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 909 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 910 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 911 912 IF(ll_tmp1) THEN 913 zcpx(ji,jj) = 1.0_wp 914 ELSE IF(ll_tmp2) THEN 915 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 916 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 917 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 918 919 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 920 ELSE 921 zcpx(ji,jj) = 0._wp 922 END IF 923 924 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 925 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 926 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 927 & > rn_wdmin1 + rn_wdmin2 928 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 929 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 930 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 931 932 IF(ll_tmp1) THEN 933 zcpy(ji,jj) = 1.0_wp 934 ELSE IF(ll_tmp2) THEN 935 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 936 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 937 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 938 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 939 940 ELSE 941 zcpy(ji,jj) = 0._wp 942 ENDIF 943 END_2D 1019 944 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. ) 1020 945 ENDIF … … 1025 950 1026 951 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 1027 DO jj = 1, jpj 1028 DO ji = 1, jpi 1029 jk = mbkt(ji,jj)+1 1030 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1031 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 1032 ELSEIF( jk < jpkm1 ) THEN 1033 DO jkk = jk+1, jpk 1034 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), & 1035 & gde3w_n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1036 END DO 1037 ENDIF 1038 END DO 1039 END DO 952 DO_2D_11_11 953 jk = mbkt(ji,jj)+1 954 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 955 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 956 ELSEIF( jk < jpkm1 ) THEN 957 DO jkk = jk+1, jpk 958 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 959 & gde3w(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 960 END DO 961 ENDIF 962 END_2D 1040 963 1041 964 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 1042 DO jj = 1, jpj 1043 DO ji = 1, jpi 1044 zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 1045 END DO 1046 END DO 1047 1048 DO jk = 2, jpk 1049 DO jj = 1, jpj 1050 DO ji = 1, jpi 1051 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 1052 END DO 1053 END DO 1054 END DO 965 DO_2D_11_11 966 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 967 END_2D 968 969 DO_3D_11_11( 2, jpk ) 970 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 971 END_3D 1055 972 1056 973 fsp(:,:,:) = zrhh (:,:,:) … … 1063 980 1064 981 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1065 DO jj = 2, jpj 1066 DO ji = 2, jpi 1067 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1068 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1069 1070 ! assuming linear profile across the top half surface layer 1071 zhpi(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) * zrhdt1 1072 END DO 1073 END DO 982 DO_2D_01_01 983 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 984 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 985 986 ! assuming linear profile across the top half surface layer 987 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 988 END_2D 1074 989 1075 990 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 1076 DO jk = 2, jpkm1 1077 DO jj = 2, jpj 1078 DO ji = 2, jpi 1079 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1080 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1081 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1082 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1083 END DO 1084 END DO 1085 END DO 991 DO_3D_01_01( 2, jpkm1 ) 992 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 993 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 994 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 995 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 996 END_3D 1086 997 1087 998 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1088 999 1089 1000 ! Prepare zsshu_n and zsshv_n 1090 DO jj = 2, jpjm1 1091 DO ji = 2, jpim1 1001 DO_2D_00_00 1092 1002 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1093 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &1003 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1094 1004 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1095 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &1005 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1096 1006 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1097 1007 !!gm not this: 1098 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1099 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1100 zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 1101 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1102 END DO 1103 END DO 1008 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1009 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1010 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1011 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1012 END_2D 1104 1013 1105 1014 CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1., zsshv_n, 'V', 1. ) 1106 1015 1107 DO jj = 2, jpjm1 1108 DO ji = 2, jpim1 1109 zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 1110 zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 1111 END DO 1112 END DO 1113 1114 DO jk = 2, jpkm1 1115 DO jj = 2, jpjm1 1116 DO ji = 2, jpim1 1117 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 1118 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1119 END DO 1120 END DO 1121 END DO 1122 1123 DO jk = 1, jpkm1 1124 DO jj = 2, jpjm1 1125 DO ji = 2, jpim1 1126 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 1127 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 1128 END DO 1129 END DO 1130 END DO 1131 1132 DO jk = 1, jpkm1 1133 DO jj = 2, jpjm1 1134 DO ji = 2, jpim1 1135 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1136 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1137 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1138 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1139 END DO 1140 END DO 1141 END DO 1142 1143 1144 DO jk = 1, jpkm1 1145 DO jj = 2, jpjm1 1146 DO ji = 2, jpim1 1147 zpwes = 0._wp; zpwed = 0._wp 1148 zpnss = 0._wp; zpnsd = 0._wp 1149 zuijk = zu(ji,jj,jk) 1150 zvijk = zv(ji,jj,jk) 1151 1152 !!!!! for u equation 1153 IF( jk <= mbku(ji,jj) ) THEN 1154 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1155 jis = ji + 1; jid = ji 1156 ELSE 1157 jis = ji; jid = ji +1 1158 ENDIF 1159 1160 ! integrate the pressure on the shallow side 1161 jk1 = jk 1162 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1163 IF( jk1 == mbku(ji,jj) ) THEN 1164 zuijk = -zdept(jis,jj,jk1) 1165 EXIT 1166 ENDIF 1167 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1168 zpwes = zpwes + & 1169 integ_spline(zdept(jis,jj,jk1), zdeps, & 1170 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1171 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1172 jk1 = jk1 + 1 1173 END DO 1174 1175 ! integrate the pressure on the deep side 1176 jk1 = jk 1177 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1178 IF( jk1 == 1 ) THEN 1179 zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 1180 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1181 bsp(jid,jj,1), csp(jid,jj,1), & 1182 dsp(jid,jj,1)) * zdeps 1183 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1184 EXIT 1185 ENDIF 1186 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1187 zpwed = zpwed + & 1188 integ_spline(zdeps, zdept(jid,jj,jk1), & 1189 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1190 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1191 jk1 = jk1 - 1 1192 END DO 1193 1194 ! update the momentum trends in u direction 1195 1196 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1197 IF( .NOT.ln_linssh ) THEN 1198 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1199 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1200 ELSE 1201 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1202 ENDIF 1203 IF( ln_wd_il ) THEN 1204 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1205 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1206 ENDIF 1207 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1208 ENDIF 1209 1210 !!!!! for v equation 1211 IF( jk <= mbkv(ji,jj) ) THEN 1212 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 1213 jjs = jj + 1; jjd = jj 1214 ELSE 1215 jjs = jj ; jjd = jj + 1 1216 ENDIF 1217 1218 ! integrate the pressure on the shallow side 1219 jk1 = jk 1220 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 1221 IF( jk1 == mbkv(ji,jj) ) THEN 1222 zvijk = -zdept(ji,jjs,jk1) 1223 EXIT 1224 ENDIF 1225 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 1226 zpnss = zpnss + & 1227 integ_spline(zdept(ji,jjs,jk1), zdeps, & 1228 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1229 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1230 jk1 = jk1 + 1 1231 END DO 1232 1233 ! integrate the pressure on the deep side 1234 jk1 = jk 1235 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1236 IF( jk1 == 1 ) THEN 1237 zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 1238 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1239 bsp(ji,jjd,1), csp(ji,jjd,1), & 1240 dsp(ji,jjd,1) ) * zdeps 1241 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1242 EXIT 1243 ENDIF 1244 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 1245 zpnsd = zpnsd + & 1246 integ_spline(zdeps, zdept(ji,jjd,jk1), & 1247 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1248 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1249 jk1 = jk1 - 1 1250 END DO 1251 1252 1253 ! update the momentum trends in v direction 1254 1255 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1256 IF( .NOT.ln_linssh ) THEN 1257 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1258 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1259 ELSE 1260 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1261 ENDIF 1262 IF( ln_wd_il ) THEN 1263 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1264 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1265 ENDIF 1266 1267 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1268 ENDIF 1269 ! 1270 END DO 1016 DO_2D_00_00 1017 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1018 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1019 END_2D 1020 1021 DO_3D_00_00( 2, jpkm1 ) 1022 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1023 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1024 END_3D 1025 1026 DO_3D_00_00( 1, jpkm1 ) 1027 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1028 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1029 END_3D 1030 1031 DO_3D_00_00( 1, jpkm1 ) 1032 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1033 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1034 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1035 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1036 END_3D 1037 1038 1039 DO_3D_00_00( 1, jpkm1 ) 1040 zpwes = 0._wp; zpwed = 0._wp 1041 zpnss = 0._wp; zpnsd = 0._wp 1042 zuijk = zu(ji,jj,jk) 1043 zvijk = zv(ji,jj,jk) 1044 1045 !!!!! for u equation 1046 IF( jk <= mbku(ji,jj) ) THEN 1047 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1048 jis = ji + 1; jid = ji 1049 ELSE 1050 jis = ji; jid = ji +1 1051 ENDIF 1052 1053 ! integrate the pressure on the shallow side 1054 jk1 = jk 1055 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1056 IF( jk1 == mbku(ji,jj) ) THEN 1057 zuijk = -zdept(jis,jj,jk1) 1058 EXIT 1059 ENDIF 1060 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1061 zpwes = zpwes + & 1062 integ_spline(zdept(jis,jj,jk1), zdeps, & 1063 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1064 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1065 jk1 = jk1 + 1 1271 1066 END DO 1272 END DO 1067 1068 ! integrate the pressure on the deep side 1069 jk1 = jk 1070 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1071 IF( jk1 == 1 ) THEN 1072 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1073 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1074 bsp(jid,jj,1), csp(jid,jj,1), & 1075 dsp(jid,jj,1)) * zdeps 1076 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1077 EXIT 1078 ENDIF 1079 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1080 zpwed = zpwed + & 1081 integ_spline(zdeps, zdept(jid,jj,jk1), & 1082 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1083 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1084 jk1 = jk1 - 1 1085 END DO 1086 1087 ! update the momentum trends in u direction 1088 1089 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1090 IF( .NOT.ln_linssh ) THEN 1091 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1092 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1093 ELSE 1094 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1095 ENDIF 1096 IF( ln_wd_il ) THEN 1097 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1098 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1099 ENDIF 1100 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1101 ENDIF 1102 1103 !!!!! for v equation 1104 IF( jk <= mbkv(ji,jj) ) THEN 1105 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 1106 jjs = jj + 1; jjd = jj 1107 ELSE 1108 jjs = jj ; jjd = jj + 1 1109 ENDIF 1110 1111 ! integrate the pressure on the shallow side 1112 jk1 = jk 1113 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 1114 IF( jk1 == mbkv(ji,jj) ) THEN 1115 zvijk = -zdept(ji,jjs,jk1) 1116 EXIT 1117 ENDIF 1118 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 1119 zpnss = zpnss + & 1120 integ_spline(zdept(ji,jjs,jk1), zdeps, & 1121 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1122 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1123 jk1 = jk1 + 1 1124 END DO 1125 1126 ! integrate the pressure on the deep side 1127 jk1 = jk 1128 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1129 IF( jk1 == 1 ) THEN 1130 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1131 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1132 bsp(ji,jjd,1), csp(ji,jjd,1), & 1133 dsp(ji,jjd,1) ) * zdeps 1134 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1135 EXIT 1136 ENDIF 1137 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 1138 zpnsd = zpnsd + & 1139 integ_spline(zdeps, zdept(ji,jjd,jk1), & 1140 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1141 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1142 jk1 = jk1 - 1 1143 END DO 1144 1145 1146 ! update the momentum trends in v direction 1147 1148 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1149 IF( .NOT.ln_linssh ) THEN 1150 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1151 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1152 ELSE 1153 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1154 ENDIF 1155 IF( ln_wd_il ) THEN 1156 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1157 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1158 ENDIF 1159 1160 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1161 ENDIF 1162 ! 1163 END_3D 1273 1164 ! 1274 1165 IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
Note: See TracChangeset
for help on using the changeset viewer.
