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