Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (53 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5930 r6140 45 45 USE wrk_nemo ! Memory Allocation 46 46 USE timing ! Timing 47 USE iom 47 48 48 49 IMPLICIT NONE … … 52 53 PUBLIC dyn_hpg_init ! routine called by opa module 53 54 54 ! 55 LOGICAL , PUBLIC :: ln_hpg_zco 56 LOGICAL , PUBLIC :: ln_hpg_zps 57 LOGICAL , PUBLIC :: ln_hpg_sco 58 LOGICAL , PUBLIC :: ln_hpg_djc 59 LOGICAL , PUBLIC :: ln_hpg_prj 60 LOGICAL , PUBLIC :: ln_hpg_isf 55 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 56 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 57 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 58 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 59 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 60 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 61 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 61 62 62 63 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 63 64 64 65 !! * Substitutions 65 # include "domzgr_substitute.h90"66 66 # include "vectopt_loop_substitute.h90" 67 67 !!---------------------------------------------------------------------- … … 131 131 INTEGER :: ios ! Local integer output status for namelist read 132 132 !! 133 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 134 REAL(wp) :: znad 135 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 136 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 137 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 138 !! 133 139 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 140 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf … … 137 143 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 138 144 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 139 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )140 145 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp ) 146 ! 141 147 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 142 148 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 143 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )149 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 144 150 IF(lwm) WRITE ( numond, namdyn_hpg ) 145 151 ! … … 162 168 & either ln_hpg_sco or ln_hpg_prj instead') 163 169 ! 164 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )&165 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&166 & the standard jacobian formulation hpg_sco or&167 & the pressure jacobian formulation hpg_prj')170 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 171 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 172 & ' the standard jacobian formulation hpg_sco or ' , & 173 & ' the pressure jacobian formulation hpg_prj' ) 168 174 169 175 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & … … 190 196 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 191 197 ! 192 ! initialisation of ice load 193 riceload(:,:)=0.0 198 ! initialisation of ice shelf load 199 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 200 IF ( ln_isfcav ) THEN 201 CALL wrk_alloc( jpi,jpj, 2, ztstop) 202 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 203 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 204 ! 205 IF(lwp) WRITE(numout,*) 206 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 207 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 208 209 ! To use density and not density anomaly 210 znad=1._wp 211 212 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 213 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 214 215 ! compute density of the water displaced by the ice shelf 216 DO jk = 1, jpk 217 CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk)) 218 END DO 219 220 ! compute rhd at the ice/oce interface (ice shelf side) 221 CALL eos(ztstop,risfdep,zrhdtop_isf) 222 223 ! Surface value + ice shelf gradient 224 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 225 ! divided by 2 later 226 ziceload = 0._wp 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ikt=mikt(ji,jj) 230 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 231 DO jk=2,ikt-1 232 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 233 & * (1._wp - tmask(ji,jj,jk)) 234 END DO 235 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 236 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 237 END DO 238 END DO 239 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 240 241 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 242 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 243 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 244 END IF 194 245 ! 195 246 END SUBROUTINE dyn_hpg_init … … 213 264 !!---------------------------------------------------------------------- 214 265 INTEGER, INTENT(in) :: kt ! ocean time-step index 215 ! !266 ! 216 267 INTEGER :: ji, jj, jk ! dummy loop indices 217 268 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars … … 219 270 !!---------------------------------------------------------------------- 220 271 ! 221 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )272 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 222 273 ! 223 274 IF( kt == nit000 ) THEN … … 232 283 DO jj = 2, jpjm1 233 284 DO ji = fs_2, fs_jpim1 ! vector opt. 234 zcoef1 = zcoef0 * fse3w(ji,jj,1)285 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 235 286 ! hydrostatic pressure gradient 236 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) /e1u(ji,jj)237 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)287 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 288 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 238 289 ! add to the general momentum trend 239 290 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 247 298 DO jj = 2, jpjm1 248 299 DO ji = fs_2, fs_jpim1 ! vector opt. 249 zcoef1 = zcoef0 * fse3w(ji,jj,jk)300 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 250 301 ! hydrostatic pressure gradient 251 302 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 252 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &253 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)303 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 304 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 254 305 255 306 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 256 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &257 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)307 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 308 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 258 309 ! add to the general momentum trend 259 310 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 263 314 END DO 264 315 ! 265 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )316 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 266 317 ! 267 318 END SUBROUTINE hpg_zco … … 284 335 !!---------------------------------------------------------------------- 285 336 ! 286 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )337 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 287 338 ! 288 339 IF( kt == nit000 ) THEN … … 301 352 DO jj = 2, jpjm1 302 353 DO ji = fs_2, fs_jpim1 ! vector opt. 303 zcoef1 = zcoef0 * fse3w(ji,jj,1)354 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 304 355 ! hydrostatic pressure gradient 305 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) /e1u(ji,jj)306 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)356 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 357 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 307 358 ! add to the general momentum trend 308 359 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 310 361 END DO 311 362 END DO 312 313 363 314 364 ! interior value (2=<jk=<jpkm1) … … 316 366 DO jj = 2, jpjm1 317 367 DO ji = fs_2, fs_jpim1 ! vector opt. 318 zcoef1 = zcoef0 * fse3w(ji,jj,jk)368 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 319 369 ! hydrostatic pressure gradient 320 370 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 321 371 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 322 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)372 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 323 373 324 374 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 325 375 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 326 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)376 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 327 377 ! add to the general momentum trend 328 378 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 331 381 END DO 332 382 END DO 333 334 383 335 384 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 338 387 iku = mbku(ji,jj) 339 388 ikv = mbkv(ji,jj) 340 zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )341 zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )389 zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj ,iku) ) 390 zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) ) 342 391 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 343 392 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 344 393 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 345 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) /e1u(ji,jj)394 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 346 395 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 347 396 ENDIF … … 349 398 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 350 399 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 351 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) /e2v(ji,jj)400 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 352 401 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 353 402 ENDIF … … 355 404 END DO 356 405 ! 357 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )406 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 358 407 ! 359 408 END SUBROUTINE hpg_zps 409 360 410 361 411 SUBROUTINE hpg_sco( kt ) … … 391 441 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 392 442 ENDIF 393 394 ! Local constant initialization 443 ! 395 444 zcoef0 = - grav * 0.5_wp 396 ! To use density and not density anomaly 397 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 398 ELSE ; znad = 0._wp ! Fixed volume 445 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 446 ELSE ; znad = 1._wp ! Variable volume: density 399 447 ENDIF 400 448 ! 401 449 ! Surface value 402 450 DO jj = 2, jpjm1 403 451 DO ji = fs_2, fs_jpim1 ! vector opt. 404 452 ! hydrostatic pressure gradient along s-surfaces 405 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) )&406 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))407 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) )&408 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))453 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 454 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 455 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 456 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 409 457 ! s-coordinate pressure gradient correction 410 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd(ji,jj,1) + 2._wp * znad ) &411 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)412 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd(ji,jj,1) + 2._wp * znad ) &413 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)458 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 459 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 460 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 461 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 414 462 ! add to the general momentum trend 415 463 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 423 471 DO ji = fs_2, fs_jpim1 ! vector opt. 424 472 ! hydrostatic pressure gradient along s-surfaces 425 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 /e1u(ji,jj) &426 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &427 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )428 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 /e2v(ji,jj) &429 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &430 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )473 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 474 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 475 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 476 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 477 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 478 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 431 479 ! s-coordinate pressure gradient correction 432 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd(ji,jj,jk) + 2._wp * znad ) &433 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) /e1u(ji,jj)434 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd(ji,jj,jk) + 2._wp * znad ) &435 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) /e2v(ji,jj)480 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 481 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 482 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 483 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 436 484 ! add to the general momentum trend 437 485 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 445 493 END SUBROUTINE hpg_sco 446 494 495 447 496 SUBROUTINE hpg_isf( kt ) 448 497 !!--------------------------------------------------------------------- 449 !! *** ROUTINE hpg_ sco***498 !! *** ROUTINE hpg_isf *** 450 499 !! 451 500 !! ** Method : s-coordinate case. Jacobian scheme. … … 466 515 INTEGER, INTENT(in) :: kt ! ocean time-step index 467 516 !! 468 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices469 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars470 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd517 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 518 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 519 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 471 520 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 472 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 473 !!---------------------------------------------------------------------- 474 ! 475 CALL wrk_alloc( jpi,jpj, 2, ztstop) 476 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 477 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 478 ! 479 IF( kt == nit000 ) THEN 480 IF(lwp) WRITE(numout,*) 481 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 482 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 483 ENDIF 484 521 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 522 !!---------------------------------------------------------------------- 523 ! 524 CALL wrk_alloc( jpi,jpj, 2, ztstop) 525 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 526 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 527 ! 485 528 ! Local constant initialization 486 529 zcoef0 = - grav * 0.5_wp 530 487 531 ! To use density and not density anomaly 488 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume489 ! ELSE ; znad = 0._wp ! Fixed volume490 ! ENDIF491 532 znad=1._wp 533 492 534 ! iniitialised to 0. zhpi zhpi 493 535 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 536 495 ! Partial steps: top & bottom before horizontal gradient496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, &497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi )499 500 !==================================================================================501 !=====Compute iceload and contribution of the half first wet layer =================502 !===================================================================================503 504 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)505 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp506 507 ! compute density of the water displaced by the ice shelf508 zrhd = rhd ! save rhd509 DO jk = 1, jpk510 zdept(:,:)=gdept_1d(jk)511 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))512 END DO513 WHERE ( tmask(:,:,:) == 1._wp)514 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd515 END WHERE516 517 ! compute rhd at the ice/oce interface (ice shelf side)518 CALL eos(ztstop,risfdep,zrhdtop_isf)519 520 537 ! compute rhd at the ice/oce interface (ocean side) 538 ! usefull to reduce residual current in the test case ISOMIP with no melting 521 539 DO ji=1,jpi 522 540 DO jj=1,jpj … … 526 544 END DO 527 545 END DO 528 CALL eos(ztstop,risfdep,zrhdtop_oce) 529 ! 530 ! Surface value + ice shelf gradient 531 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 532 ziceload = 0._wp 533 DO jj = 1, jpj 534 DO ji = 1, jpi ! vector opt. 535 ikt=mikt(ji,jj) 536 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 537 DO jk=2,ikt-1 538 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 539 & * (1._wp - tmask(ji,jj,jk)) 540 END DO 541 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 542 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 543 END DO 544 END DO 545 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 546 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 546 CALL eos( ztstop, risfdep, zrhdtop_oce ) 547 548 !================================================================================== 549 !===== Compute surface value ===================================================== 550 !================================================================================== 547 551 DO jj = 2, jpjm1 548 552 DO ji = fs_2, fs_jpim1 ! vector opt. 549 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 553 ikt = mikt(ji,jj) 554 iktp1i = mikt(ji+1,jj) 555 iktp1j = mikt(ji,jj+1) 550 556 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 551 557 ! we assume ISF is in isostatic equilibrium 552 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i) &553 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj) ) &554 & - 0.5_wp * fse3w(ji ,jj ,ikt )&555 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&556 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)))557 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j) &558 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) &559 & - 0.5_wp * fse3w(ji ,jj ,ikt )&560 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&561 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ))558 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i) & 559 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 560 & - 0.5_wp * e3w_n(ji,jj,ikt) & 561 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 562 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 563 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j) & 564 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 565 & - 0.5_wp * e3w_n(ji,jj,ikt) & 566 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 567 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 562 568 ! s-coordinate pressure gradient correction (=0 if z coordinate) 563 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd(ji,jj,1) + 2._wp * znad ) &564 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)565 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd(ji,jj,1) + 2._wp * znad ) &566 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)569 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 570 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 571 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 572 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 567 573 ! add to the general momentum trend 568 574 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) … … 571 577 END DO 572 578 !================================================================================== 573 !===== Compute partial cell contribution for the top cell =========================574 !==================================================================================575 DO jj = 2, jpjm1576 DO ji = fs_2, fs_jpim1 ! vector opt.577 iku = miku(ji,jj) ;578 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp579 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))580 ! u direction581 IF ( iku .GT. 1 ) THEN582 ! case iku583 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &584 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &585 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )586 ! corrective term ( = 0 if z coordinate )587 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)588 ! zhpi will be added in interior loop589 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap590 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure591 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)592 593 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)594 zhpiint = zcoef0 / e1u(ji,jj) &595 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &596 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &597 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &598 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )599 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint600 END IF601 602 ! v direction603 ikv = mikv(ji,jj)604 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))605 IF ( ikv .GT. 1 ) THEN606 ! case ikv607 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &608 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &609 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )610 ! corrective term ( = 0 if z coordinate )611 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)612 ! zhpi will be added in interior loop613 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap614 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure615 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)616 617 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)618 zhpjint = zcoef0 / e2v(ji,jj) &619 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &620 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &621 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &622 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )623 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint624 END IF625 END DO626 END DO627 628 !==================================================================================629 579 !===== Compute interior value ===================================================== 630 580 !================================================================================== 631 632 DO jj = 2, jpjm1 633 DO ji = fs_2, fs_jpim1 ! vector opt. 634 iku=miku(ji,jj); ikv=mikv(ji,jj) 635 DO jk = 2, jpkm1 581 ! interior value (2=<jk=<jpkm1) 582 DO jk = 2, jpkm1 583 DO jj = 2, jpjm1 584 DO ji = fs_2, fs_jpim1 ! vector opt. 636 585 ! hydrostatic pressure gradient along s-surfaces 637 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 638 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 639 & + zcoef0 / e1u(ji,jj) & 640 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 641 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 642 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 643 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 586 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 587 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 588 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 589 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 590 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 591 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 644 592 ! s-coordinate pressure gradient correction 645 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 646 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 647 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 648 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 649 650 ! hydrostatic pressure gradient along s-surfaces 651 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 652 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 653 & + zcoef0 / e2v(ji,jj) & 654 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 655 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 656 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 657 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 658 ! s-coordinate pressure gradient correction 659 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 660 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 661 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 593 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 594 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 595 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 596 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 662 597 ! add to the general momentum trend 663 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 598 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 599 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 664 600 END DO 665 601 END DO 666 602 END DO 667 668 !================================================================================== 669 !===== Compute bottom cell contribution (partial cell) ============================ 670 !================================================================================== 671 672 DO jj = 2, jpjm1 673 DO ji = 2, jpim1 674 iku = mbku(ji,jj) 675 ikv = mbkv(ji,jj) 676 677 IF (iku .GT. 1) THEN 678 ! remove old value (interior case) 679 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 680 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 681 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 682 ! put new value 683 ! -zpshpi to avoid double contribution of the partial step in the top layer 684 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 685 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 686 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 687 END IF 688 ! v direction 689 IF (ikv .GT. 1) THEN 690 ! remove old value (interior case) 691 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 692 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 693 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 694 ! put new value 695 ! -zpshpj to avoid double contribution of the partial step in the top layer 696 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 697 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 698 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 699 END IF 700 END DO 701 END DO 702 703 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 704 rhd = zrhd 705 ! 706 CALL wrk_dealloc( jpi,jpj,2, ztstop) 707 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 708 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 603 ! 604 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 605 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 606 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 709 607 ! 710 608 END SUBROUTINE hpg_isf … … 756 654 DO jj = 2, jpjm1 757 655 DO ji = fs_2, fs_jpim1 ! vector opt. 758 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd(ji,jj,jk-1)759 dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)760 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd(ji,jj,jk )761 dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )762 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd(ji,jj,jk )763 dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )656 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 657 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 658 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 659 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 660 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 661 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 764 662 END DO 765 663 END DO … … 843 741 !------------------------------------------------------------- 844 742 845 !!bug gm : e3w- de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified846 ! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be743 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified 744 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 847 745 848 746 DO jj = 2, jpjm1 849 747 DO ji = fs_2, fs_jpim1 ! vector opt. 850 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &851 & * ( rhd(ji,jj,1) &852 & + 0.5_wp * ( rhd (ji,jj,2) - rhd(ji,jj,1) )&853 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )&854 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )748 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) & 749 & * ( rhd(ji,jj,1) & 750 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 751 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) & 752 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) ) 855 753 END DO 856 754 END DO … … 863 761 DO ji = fs_2, fs_jpim1 ! vector opt. 864 762 865 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd(ji,jj,jk-1) ) &866 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &867 & - grav * z1_10 * ( 868 & ( drhow (ji,jj,jk) - drhow(ji,jj,jk-1) ) &869 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &870 & - ( dzw (ji,jj,jk) - dzw(ji,jj,jk-1) ) &871 & * ( rhd (ji,jj,jk) - rhd(ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &763 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 764 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) & 765 & - grav * z1_10 * ( & 766 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 767 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 768 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 769 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 872 770 & ) 873 771 874 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd(ji,jj,jk) ) &875 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &876 & - grav* z1_10 * ( 877 & ( drhou (ji+1,jj,jk) - drhou(ji,jj,jk) ) &878 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &879 & - ( dzu (ji+1,jj,jk) - dzu(ji,jj,jk) ) &880 & * ( rhd (ji+1,jj,jk) - rhd(ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &772 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 773 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) & 774 & - grav* z1_10 * ( & 775 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 776 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 777 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 778 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 881 779 & ) 882 780 883 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) )&884 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &885 & - grav* z1_10 * ( 886 & ( drhov (ji,jj+1,jk) - drhov(ji,jj,jk) ) &887 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &888 & - ( dzv (ji,jj+1,jk) - dzv(ji,jj,jk) ) &889 & * ( rhd (ji,jj+1,jk) - rhd(ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &781 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 782 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) & 783 & - grav* z1_10 * ( & 784 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 785 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 786 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 787 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 890 788 & ) 891 789 … … 903 801 DO jj = 2, jpjm1 904 802 DO ji = fs_2, fs_jpim1 ! vector opt. 905 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) /e1u(ji,jj)906 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) /e2v(ji,jj)803 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 804 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 907 805 ! add to the general momentum trend 908 806 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 920 818 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 921 819 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 922 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) /e1u(ji,jj)820 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 923 821 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 924 822 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 925 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) /e2v(ji,jj)823 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 926 824 ! add to the general momentum trend 927 825 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 953 851 !! 954 852 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 955 REAL(wp) :: zcoef0, znad ! temporaryscalars956 ! !853 REAL(wp) :: zcoef0, znad ! local scalars 854 ! 957 855 !! The local variables for the correction term 958 856 INTEGER :: jk1, jis, jid, jjs, jjd … … 965 863 !!---------------------------------------------------------------------- 966 864 ! 967 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )968 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )969 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )865 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 866 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 867 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 970 868 ! 971 869 IF( kt == nit000 ) THEN … … 975 873 ENDIF 976 874 977 !!----------------------------------------------------------------------978 875 ! Local constant initialization 979 876 zcoef0 = - grav 980 znad = 0.0_wp981 IF( l k_vvl ) znad = 1._wp877 znad = 1._wp 878 IF( ln_linssh ) znad = 0._wp 982 879 983 880 ! Clean 3-D work arrays … … 989 886 DO ji = 1, jpi 990 887 jk = mbathy(ji,jj) 991 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp992 ELSE IF(jk == 1) THEN; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)993 ELSE IF(jk < jpkm1) THEN888 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 889 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 890 ELSEIF( jk < jpkm1 ) THEN 994 891 DO jkk = jk+1, jpk 995 zrhh(ji,jj,jkk) = interp1( fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),&996 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))892 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), & 893 & gde3w_n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 997 894 END DO 998 895 ENDIF … … 1003 900 DO jj = 1, jpj 1004 901 DO ji = 1, jpi 1005 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad902 zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 1006 903 END DO 1007 904 END DO … … 1010 907 DO jj = 1, jpj 1011 908 DO ji = 1, jpi 1012 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)909 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 1013 910 END DO 1014 911 END DO … … 1021 918 ! constrained cubic spline interpolation 1022 919 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 1023 CALL cspline( fsp,xsp,asp,bsp,csp,dsp,polynomial_type)920 CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1024 921 1025 922 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1026 923 DO jj = 2, jpj 1027 924 DO ji = 2, jpi 1028 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 1029 bsp(ji,jj,1), csp(ji,jj,1), & 1030 dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1) 925 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 926 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1031 927 1032 928 ! assuming linear profile across the top half surface layer 1033 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1929 zhpi(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) * zrhdt1 1034 930 END DO 1035 931 END DO … … 1039 935 DO jj = 2, jpj 1040 936 DO ji = 2, jpi 1041 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1042 integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&1043 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), &1044 csp(ji,jj,jk-1), dsp(ji,jj,jk-1))937 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 938 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 939 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 940 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1045 941 END DO 1046 942 END DO … … 1052 948 DO jj = 2, jpjm1 1053 949 DO ji = 2, jpim1 950 !!gm BUG ? if it is ssh at u- & v-point then it should be: 951 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 952 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 953 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 954 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 955 !!gm not this: 1054 956 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1055 957 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp … … 1061 963 DO jj = 2, jpjm1 1062 964 DO ji = 2, jpim1 1063 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)1064 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad)965 zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 966 zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 1065 967 END DO 1066 968 END DO … … 1069 971 DO jj = 2, jpjm1 1070 972 DO ji = 2, jpim1 1071 zu(ji,jj,jk) = zu(ji,jj,jk-1) - fse3u(ji,jj,jk)1072 zv(ji,jj,jk) = zv(ji,jj,jk-1) - fse3v(ji,jj,jk)973 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 974 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1073 975 END DO 1074 976 END DO … … 1078 980 DO jj = 2, jpjm1 1079 981 DO ji = 2, jpim1 1080 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)1081 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)982 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 983 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 1082 984 END DO 1083 985 END DO … … 1087 989 DO jj = 2, jpjm1 1088 990 DO ji = 2, jpim1 1089 zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1090 zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1091 zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1092 zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))991 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 992 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 993 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 994 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1093 995 END DO 1094 996 END DO … … 1148 1050 ! update the momentum trends in u direction 1149 1051 1150 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))1151 IF( lk_vvl) THEN1152 zdpdx2 = zcoef0 /e1u(ji,jj) * &1153 1052 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1053 IF( .NOT.ln_linssh ) THEN 1054 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1055 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1154 1056 ELSE 1155 zdpdx2 = zcoef0 /e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)1057 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1156 1058 ENDIF 1157 1158 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1159 &umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)1059 !!gm Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:) by definition 1060 !!gm in the line below only * umask(ji,jj,jk) is needed !! 1061 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1160 1062 ENDIF 1161 1063 … … 1205 1107 ! update the momentum trends in v direction 1206 1108 1207 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))1208 IF( lk_vvl) THEN1209 zdpdy2 = zcoef0 /e2v(ji,jj) * &1210 1109 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1110 IF( .NOT.ln_linssh ) THEN 1111 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1112 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1211 1113 ELSE 1212 zdpdy2 = zcoef0 /e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1114 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1213 1115 ENDIF 1214 1215 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1216 &vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)1116 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition 1117 !!gm in the line below only * vmask(ji,jj,jk) is needed !! 1118 va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1217 1119 ENDIF 1218 1219 1220 END DO 1221 END DO 1222 END DO 1223 ! 1224 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1225 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1226 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1120 ! 1121 END DO 1122 END DO 1123 END DO 1124 ! 1125 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1126 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1127 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1227 1128 ! 1228 1129 END SUBROUTINE hpg_prj 1229 1130 1230 1131 1231 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type)1132 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1232 1133 !!---------------------------------------------------------------------- 1233 1134 !! *** ROUTINE cspline *** … … 1239 1140 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1240 1141 !!---------------------------------------------------------------------- 1241 IMPLICIT NONE 1242 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1243 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1244 ! the interpoated function 1245 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1246 ! 2: Linear 1142 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate 1143 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1144 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1247 1145 ! 1248 1146 INTEGER :: ji, jj, jk ! dummy loop indices … … 1252 1150 REAL(wp) :: zdf(size(fsp,3)) 1253 1151 !!---------------------------------------------------------------------- 1254 1152 ! 1153 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1255 1154 jpi = size(fsp,1) 1256 1155 jpj = size(fsp,2) 1257 1156 jpkm1 = size(fsp,3) - 1 1258 1259 1157 ! 1260 1158 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1261 1159 DO ji = 1, jpi … … 1290 1188 1291 1189 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1292 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)1190 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1293 1191 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1294 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1295 & 0.5_wp * zdf(jpkm1 - 1) 1192 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1296 1193 1297 1194 DO jk = 1, jpkm1 - 1 … … 1316 1213 END DO 1317 1214 1318 ELSE IF (polynomial_type == 2) THEN ! Linear1215 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1319 1216 DO ji = 1, jpi 1320 1217 DO jj = 1, jpj … … 1347 1244 !! extrapolation is also permitted (no value limit) 1348 1245 !!---------------------------------------------------------------------- 1349 IMPLICIT NONE1350 1246 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1351 1247 REAL(wp) :: f ! result of the interpolation (extrapolation) 1352 1248 REAL(wp) :: zdeltx 1353 1249 !!---------------------------------------------------------------------- 1354 1250 ! 1355 1251 zdeltx = xr - xl 1356 IF( abs(zdeltx) <= 10._wp * EPSILON(x)) THEN1357 f = 0.5_wp * (fl + fr)1252 IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 1253 f = 0.5_wp * (fl + fr) 1358 1254 ELSE 1359 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx1255 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1360 1256 ENDIF 1361 1257 ! 1362 1258 END FUNCTION interp1 1363 1259 1364 1260 1365 FUNCTION interp2( x, a, b, c, d) RESULT(f)1261 FUNCTION interp2( x, a, b, c, d ) RESULT(f) 1366 1262 !!---------------------------------------------------------------------- 1367 1263 !! *** ROUTINE interp1 *** … … 1372 1268 !! 1373 1269 !!---------------------------------------------------------------------- 1374 IMPLICIT NONE1375 1270 REAL(wp), INTENT(in) :: x, a, b, c, d 1376 1271 REAL(wp) :: f ! value from the interpolation 1377 1272 !!---------------------------------------------------------------------- 1378 1273 ! 1379 1274 f = a + x* ( b + x * ( c + d * x ) ) 1380 1275 ! 1381 1276 END FUNCTION interp2 1382 1277 1383 1278 1384 FUNCTION interp3( x, a, b, c, d) RESULT(f)1279 FUNCTION interp3( x, a, b, c, d ) RESULT(f) 1385 1280 !!---------------------------------------------------------------------- 1386 1281 !! *** ROUTINE interp1 *** … … 1392 1287 !! 1393 1288 !!---------------------------------------------------------------------- 1394 IMPLICIT NONE1395 1289 REAL(wp), INTENT(in) :: x, a, b, c, d 1396 1290 REAL(wp) :: f ! value from the interpolation 1397 1291 !!---------------------------------------------------------------------- 1398 1292 ! 1399 1293 f = b + x * ( 2._wp * c + 3._wp * d * x) 1400 1294 ! 1401 1295 END FUNCTION interp3 1402 1296 1403 1297 1404 FUNCTION integ_spline( xl, xr, a, b, c, d) RESULT(f)1298 FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) 1405 1299 !!---------------------------------------------------------------------- 1406 1300 !! *** ROUTINE interp1 *** … … 1411 1305 !! 1412 1306 !!---------------------------------------------------------------------- 1413 IMPLICIT NONE1414 1307 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1415 1308 REAL(wp) :: za1, za2, za3 1416 1309 REAL(wp) :: f ! integration result 1417 1310 !!---------------------------------------------------------------------- 1418 1311 ! 1419 1312 za1 = 0.5_wp * b 1420 1313 za2 = c / 3.0_wp 1421 1314 za3 = 0.25_wp * d 1422 1315 ! 1423 1316 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1424 1317 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1425 1318 ! 1426 1319 END FUNCTION integ_spline 1427 1320
Note: See TracChangeset
for help on using the changeset viewer.
