- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5930 r6060 52 52 PUBLIC dyn_hpg_init ! routine called by opa module 53 53 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 54 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 55 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 56 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 57 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 58 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 59 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 60 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 61 61 62 62 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 63 63 64 64 !! * Substitutions 65 # include "domzgr_substitute.h90"66 65 # include "vectopt_loop_substitute.h90" 67 66 !!---------------------------------------------------------------------- … … 137 136 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 138 137 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 138 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp ) 139 ! 141 140 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 142 141 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 143 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )142 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 144 143 IF(lwm) WRITE ( numond, namdyn_hpg ) 145 144 ! … … 162 161 & either ln_hpg_sco or ln_hpg_prj instead') 163 162 ! 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')163 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 164 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 165 & ' the standard jacobian formulation hpg_sco or ' , & 166 & ' the pressure jacobian formulation hpg_prj' ) 168 167 169 168 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & … … 213 212 !!---------------------------------------------------------------------- 214 213 INTEGER, INTENT(in) :: kt ! ocean time-step index 215 ! !214 ! 216 215 INTEGER :: ji, jj, jk ! dummy loop indices 217 216 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars … … 219 218 !!---------------------------------------------------------------------- 220 219 ! 221 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )220 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 222 221 ! 223 222 IF( kt == nit000 ) THEN … … 232 231 DO jj = 2, jpjm1 233 232 DO ji = fs_2, fs_jpim1 ! vector opt. 234 zcoef1 = zcoef0 * fse3w(ji,jj,1)233 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 235 234 ! 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)235 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 236 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 238 237 ! add to the general momentum trend 239 238 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 247 246 DO jj = 2, jpjm1 248 247 DO ji = fs_2, fs_jpim1 ! vector opt. 249 zcoef1 = zcoef0 * fse3w(ji,jj,jk)248 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 250 249 ! hydrostatic pressure gradient 251 250 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)251 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 252 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 254 253 255 254 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)255 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 256 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 258 257 ! add to the general momentum trend 259 258 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 263 262 END DO 264 263 ! 265 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )264 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 266 265 ! 267 266 END SUBROUTINE hpg_zco … … 284 283 !!---------------------------------------------------------------------- 285 284 ! 286 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )285 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 287 286 ! 288 287 IF( kt == nit000 ) THEN … … 301 300 DO jj = 2, jpjm1 302 301 DO ji = fs_2, fs_jpim1 ! vector opt. 303 zcoef1 = zcoef0 * fse3w(ji,jj,1)302 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 304 303 ! 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)304 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 305 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 307 306 ! add to the general momentum trend 308 307 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 310 309 END DO 311 310 END DO 312 313 311 314 312 ! interior value (2=<jk=<jpkm1) … … 316 314 DO jj = 2, jpjm1 317 315 DO ji = fs_2, fs_jpim1 ! vector opt. 318 zcoef1 = zcoef0 * fse3w(ji,jj,jk)316 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 319 317 ! hydrostatic pressure gradient 320 318 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 321 319 & + 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)320 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 323 321 324 322 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 325 323 & + 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)324 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 327 325 ! add to the general momentum trend 328 326 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 331 329 END DO 332 330 END DO 333 334 331 335 332 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 338 335 iku = mbku(ji,jj) 339 336 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) )337 zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj ,iku) ) 338 zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) ) 342 339 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 343 340 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 344 341 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)342 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 346 343 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 347 344 ENDIF … … 349 346 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 350 347 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)348 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 352 349 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 353 350 ENDIF … … 355 352 END DO 356 353 ! 357 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )354 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 358 355 ! 359 356 END SUBROUTINE hpg_zps 357 360 358 361 359 SUBROUTINE hpg_sco( kt ) … … 391 389 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 392 390 ENDIF 393 394 ! Local constant initialization 391 ! 395 392 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 393 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 394 ELSE ; znad = 1._wp ! Variable volume: density 399 395 ENDIF 400 396 ! 401 397 ! Surface value 402 398 DO jj = 2, jpjm1 403 399 DO ji = fs_2, fs_jpim1 ! vector opt. 404 400 ! 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) ))401 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 402 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 403 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 404 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 409 405 ! 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)406 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 407 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 408 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 409 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 414 410 ! add to the general momentum trend 415 411 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 423 419 DO ji = fs_2, fs_jpim1 ! vector opt. 424 420 ! 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 ) )421 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 422 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 423 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 424 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 425 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 426 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 431 427 ! 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)428 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 429 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 430 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 431 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 436 432 ! add to the general momentum trend 437 433 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 444 440 ! 445 441 END SUBROUTINE hpg_sco 442 446 443 447 444 SUBROUTINE hpg_isf( kt ) … … 473 470 !!---------------------------------------------------------------------- 474 471 ! 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 ) THEN472 CALL wrk_alloc( jpi,jpj, 2, ztstop ) 473 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 474 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj ) 475 ! 476 IF( kt == nit000 ) THEN 480 477 IF(lwp) WRITE(numout,*) 481 478 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 482 479 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 483 480 ENDIF 484 485 ! Local constant initialization 481 ! 486 482 zcoef0 = - grav * 0.5_wp 487 ! To use density and not density anomaly 488 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 489 ! ELSE ; znad = 0._wp ! Fixed volume 490 ! ENDIF 491 znad=1._wp 492 ! iniitialised to 0. zhpi zhpi 493 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 495 ! Partial steps: top & bottom before horizontal gradient 496 !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 ) 483 IF( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 484 ELSE ; znad = 1._wp ! Variable volume: density 485 ENDIF 486 zhpi(:,:,:) = 0._wp 487 zhpj(:,:,:) = 0._wp 499 488 500 489 !================================================================================== … … 503 492 504 493 ! 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_wp 494 ztstop(:,:,jp_tem) = -1.9_wp 495 ztstop(:,:,jp_sal) = 34.4_wp 496 497 !!gm I have the feeling that a much simplier and faster computation can be performed... 498 !!gm ====>>>> We have to discuss ! 499 500 !!gm below, faster to compute the ISF density in zrhd and remplace rhd value where tmask=0 501 !!gm furthermore, this calculation does not depends on time : do it at the first time-step only.... 506 502 507 503 ! compute density of the water displaced by the ice shelf 508 zrhd = rhd! save rhd504 zrhd(:,:,:) = rhd(:,:,:) ! save rhd 509 505 DO jk = 1, jpk 510 zdept(:,:)=gdept_1d(jk)511 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))512 END DO 513 WHERE ( tmask(:,:,:) == 1._wp)506 zdept(:,:) = gdept_1d(jk) 507 CALL eos( ztstop(:,:,:), zdept(:,:), rhd(:,:,jk) ) 508 END DO 509 WHERE( tmask(:,:,:) == 1._wp ) 514 510 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 515 511 END WHERE 516 512 517 513 ! compute rhd at the ice/oce interface (ice shelf side) 518 CALL eos( ztstop,risfdep,zrhdtop_isf)514 CALL eos( ztstop, risfdep, zrhdtop_isf ) 519 515 520 516 ! compute rhd at the ice/oce interface (ocean side) 521 DO ji =1,jpi522 DO jj =1,jpj523 ikt =mikt(ji,jj)524 ztstop(ji,jj, 1)=tsn(ji,jj,ikt,1)525 ztstop(ji,jj, 2)=tsn(ji,jj,ikt,2)517 DO ji = 1, jpi 518 DO jj = 1, jpj 519 ikt = mikt(ji,jj) 520 ztstop(ji,jj,jp_tem) = tsn(ji,jj,ikt,jp_tem) 521 ztstop(ji,jj,jp_sal) = tsn(ji,jj,ikt,jp_sal) 526 522 END DO 527 523 END DO 528 CALL eos( ztstop,risfdep,zrhdtop_oce)524 CALL eos( ztstop, risfdep, zrhdtop_oce ) 529 525 ! 530 526 ! Surface value + ice shelf gradient … … 533 529 DO jj = 1, jpj 534 530 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-1538 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) &531 ikt = mikt(ji,jj) 532 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 533 DO jk = 2, ikt-1 534 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 539 535 & * (1._wp - tmask(ji,jj,jk)) 540 536 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 diaar5537 IF( ikt >= 2 ) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 538 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 539 END DO 540 END DO 541 riceload(:,:) = 0._wp ; riceload(:,:) = ziceload(:,:) ! need to be saved for diaar5 546 542 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 547 543 DO jj = 2, jpjm1 548 544 DO ji = fs_2, fs_jpim1 ! vector opt. 549 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 545 ikt = mikt(ji,jj) 546 iktp1i = mikt(ji+1,jj) 547 iktp1j = mikt(ji,jj+1) 550 548 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 551 549 ! 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) ) ) 550 zhpi(ji,jj,1) = zcoef0 * ( & 551 & 0.5_wp * e3w_n(ji+1,jj,iktp1i) * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 552 & - 0.5_wp * e3w_n(ji ,jj,ikt ) * ( 2._wp * znad + rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) & 553 & + ( ziceload(ji+1,jj) - ziceload(ji,jj) ) ) * r1_e1u(ji,jj) 554 zhpj(ji,jj,1) = zcoef0 * ( & 555 & 0.5_wp * e3w_n(ji,jj+1,iktp1j) * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 556 & - 0.5_wp * e3w_n(ji,jj ,ikt ) * ( 2._wp * znad + rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) & 557 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) * r1_e2v(ji,jj) 562 558 ! 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)559 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 560 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 561 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 562 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 567 563 ! add to the general momentum trend 568 564 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) … … 575 571 DO jj = 2, jpjm1 576 572 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_wp 573 iku = miku(ji,jj) 574 zpshpi(ji,jj) = 0._wp 575 zpshpj(ji,jj) = 0._wp 579 576 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 577 ! u direction 581 IF ( iku .GT.1 ) THEN578 IF( iku > 1 ) THEN 582 579 ! case iku 583 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 )580 zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu & 581 & * ( rhd(ji+1,jj,iku) + rhd(ji,jj,iku) & 582 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 586 583 ! corrective term ( = 0 if z coordinate ) 587 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) /e1u(ji,jj)584 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 588 585 ! zhpi will be added in interior loop 589 ua(ji,jj,iku) 586 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 590 587 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 591 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)= zhpi(ji,jj,iku)588 IF( mbku(ji,jj) == iku + 1 ) zpshpi(ji,jj) = zhpi(ji,jj,iku) 592 589 593 590 ! 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 & 597 & - fse3w(ji ,jj ,iku+1) * ((rhd(ji ,jj,iku+1) + znad) &598 & 599 zhpi(ji,jj,iku+1) = zcoef0 /e1u(ji,jj) * ge3rui(ji,jj) - zhpiint591 zhpiint = zcoef0 * r1_e1u(ji,jj) & 592 & * ( e3w_n(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 593 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 594 & - e3w_n(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 595 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 596 zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 600 597 END IF 601 598 … … 603 600 ikv = mikv(ji,jj) 604 601 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 ) THEN602 IF( ikv > 1 ) THEN 606 603 ! case ikv 607 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 )604 zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv & 605 & * ( rhd(ji,jj+1,ikv) + rhd(ji,jj,ikv) & 606 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 610 607 ! corrective term ( = 0 if z coordinate ) 611 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) /e2v(ji,jj)608 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 612 609 ! zhpi will be added in interior loop 613 va(ji,jj,ikv) 610 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 614 611 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 615 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) =zhpj(ji,jj,ikv)612 IF( mbkv(ji,jj) == ikv + 1 ) zpshpj(ji,jj) = zhpj(ji,jj,ikv) 616 613 617 614 ! 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 & 623 zhpj(ji,jj,ikv+1) = zcoef0 /e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint624 END 615 zhpjint = zcoef0 * r1_e2v(ji,jj) & 616 & * ( e3w_n(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 617 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 618 & - e3w_n(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 619 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 620 zhpj(ji,jj,ikv+1) = zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 621 ENDIF 625 622 END DO 626 623 END DO … … 632 629 DO jj = 2, jpjm1 633 630 DO ji = fs_2, fs_jpim1 ! vector opt. 634 iku=miku(ji,jj); ikv=mikv(ji,jj)635 631 DO jk = 2, jpkm1 636 632 ! hydrostatic pressure gradient along s-surfaces 637 633 ! 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 & 642 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad)&643 & 634 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 635 & + zcoef0 * r1_e1u(ji,jj) & 636 & * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 637 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 638 & - e3w_n(ji ,jj,jk) * ( (rhd(ji ,jj,jk ) + znad) & 639 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 644 640 ! s-coordinate pressure gradient correction 645 641 ! 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)642 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 643 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 648 644 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 649 645 650 646 ! hydrostatic pressure gradient along s-surfaces 651 647 ! 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 & 656 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad)&657 & 648 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 649 & + zcoef0 * r1_e2v(ji,jj) & 650 & * ( e3w_n(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 651 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 652 & - e3w_n(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 653 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 658 654 ! s-coordinate pressure gradient correction 659 655 ! 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)656 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 657 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 662 658 ! add to the general momentum trend 663 659 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) … … 677 673 IF (iku .GT. 1) THEN 678 674 ! 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)675 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 676 & * ( gde3w_n(ji+1,jj ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 681 677 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 682 678 ! put new value 683 679 ! -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)680 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) * r1_e1u(ji,jj) 681 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 686 682 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 687 683 END IF … … 689 685 IF (ikv .GT. 1) THEN 690 686 ! 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)687 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 688 & * ( gde3w_n(ji ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj) 693 689 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 694 690 ! put new value 695 691 ! -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)692 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj) 693 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 698 694 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 699 695 END IF … … 756 752 DO jj = 2, jpjm1 757 753 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 )754 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 755 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 756 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 757 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 758 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 759 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 764 760 END DO 765 761 END DO … … 843 839 !------------------------------------------------------------- 844 840 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 be841 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified 842 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 847 843 848 844 DO jj = 2, jpjm1 849 845 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) ) )846 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) & 847 & * ( rhd(ji,jj,1) & 848 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 849 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) & 850 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) ) 855 851 END DO 856 852 END DO … … 863 859 DO ji = fs_2, fs_jpim1 ! vector opt. 864 860 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) ) ) &861 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 862 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) & 863 & - grav * z1_10 * ( & 864 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 865 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 866 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 867 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 872 868 & ) 873 869 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) ) ) &870 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 871 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) & 872 & - grav* z1_10 * ( & 873 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 874 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 875 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 876 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 881 877 & ) 882 878 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) ) ) &879 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 880 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) & 881 & - grav* z1_10 * ( & 882 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 883 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 884 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 885 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 890 886 & ) 891 887 … … 903 899 DO jj = 2, jpjm1 904 900 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)901 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 902 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 903 ! add to the general momentum trend 908 904 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 920 916 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 921 917 & + ( ( 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)918 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 923 919 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 924 920 & + ( ( 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)921 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 926 922 ! add to the general momentum trend 927 923 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 953 949 !! 954 950 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 955 REAL(wp) :: zcoef0, znad ! temporaryscalars956 ! !951 REAL(wp) :: zcoef0, znad ! local scalars 952 ! 957 953 !! The local variables for the correction term 958 954 INTEGER :: jk1, jis, jid, jjs, jjd … … 965 961 !!---------------------------------------------------------------------- 966 962 ! 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 )963 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 964 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 965 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 970 966 ! 971 967 IF( kt == nit000 ) THEN … … 975 971 ENDIF 976 972 977 !!----------------------------------------------------------------------978 973 ! Local constant initialization 979 974 zcoef0 = - grav 980 znad = 0.0_wp981 IF( l k_vvl ) znad = 1._wp975 znad = 1._wp 976 IF( ln_linssh ) znad = 0._wp 982 977 983 978 ! Clean 3-D work arrays … … 989 984 DO ji = 1, jpi 990 985 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) THEN986 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 987 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 988 ELSEIF( jk < jpkm1 ) THEN 994 989 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))990 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), & 991 & gde3w_n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 997 992 END DO 998 993 ENDIF … … 1003 998 DO jj = 1, jpj 1004 999 DO ji = 1, jpi 1005 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad1000 zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 1006 1001 END DO 1007 1002 END DO … … 1010 1005 DO jj = 1, jpj 1011 1006 DO ji = 1, jpi 1012 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)1007 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 1013 1008 END DO 1014 1009 END DO … … 1021 1016 ! constrained cubic spline interpolation 1022 1017 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 1023 CALL cspline( fsp,xsp,asp,bsp,csp,dsp,polynomial_type)1018 CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1024 1019 1025 1020 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1026 1021 DO jj = 2, jpj 1027 1022 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) 1023 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1024 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1031 1025 1032 1026 ! assuming linear profile across the top half surface layer 1033 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt11027 zhpi(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) * zrhdt1 1034 1028 END DO 1035 1029 END DO … … 1039 1033 DO jj = 2, jpj 1040 1034 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))1035 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1036 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1037 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1038 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1045 1039 END DO 1046 1040 END DO … … 1052 1046 DO jj = 2, jpjm1 1053 1047 DO ji = 2, jpim1 1048 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1049 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 1050 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1051 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 1052 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1053 !!gm not this: 1054 1054 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1055 1055 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp … … 1061 1061 DO jj = 2, jpjm1 1062 1062 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)1063 zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 1064 zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 1065 1065 END DO 1066 1066 END DO … … 1069 1069 DO jj = 2, jpjm1 1070 1070 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)1071 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 1072 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1073 1073 END DO 1074 1074 END DO … … 1078 1078 DO jj = 2, jpjm1 1079 1079 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)1080 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 1081 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 1082 1082 END DO 1083 1083 END DO … … 1087 1087 DO jj = 2, jpjm1 1088 1088 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)))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) ) ) 1093 1093 END DO 1094 1094 END DO … … 1148 1148 ! update the momentum trends in u direction 1149 1149 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 1150 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1151 IF( .NOT.ln_linssh ) THEN 1152 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1153 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1154 1154 ELSE 1155 zdpdx2 = zcoef0 /e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)1155 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1156 1156 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)1157 !!gm Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:) by definition 1158 !!gm in the line below only * umask(ji,jj,jk) is needed !! 1159 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1160 1160 ENDIF 1161 1161 … … 1205 1205 ! update the momentum trends in v direction 1206 1206 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 1207 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1208 IF( .NOT.ln_linssh ) THEN 1209 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1210 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1211 1211 ELSE 1212 zdpdy2 = zcoef0 /e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1212 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1213 1213 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)1214 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition 1215 !!gm in the line below only * vmask(ji,jj,jk) is needed !! 1216 va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1217 1217 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 ) 1218 ! 1219 END DO 1220 END DO 1221 END DO 1222 ! 1223 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1224 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1225 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1227 1226 ! 1228 1227 END SUBROUTINE hpg_prj 1229 1228 1230 1229 1231 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type)1230 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1232 1231 !!---------------------------------------------------------------------- 1233 1232 !! *** ROUTINE cspline *** … … 1239 1238 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1240 1239 !!---------------------------------------------------------------------- 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 1240 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate 1241 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1242 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1247 1243 ! 1248 1244 INTEGER :: ji, jj, jk ! dummy loop indices … … 1252 1248 REAL(wp) :: zdf(size(fsp,3)) 1253 1249 !!---------------------------------------------------------------------- 1254 1250 ! 1251 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1255 1252 jpi = size(fsp,1) 1256 1253 jpj = size(fsp,2) 1257 1254 jpkm1 = size(fsp,3) - 1 1258 1259 1255 ! 1260 1256 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1261 1257 DO ji = 1, jpi … … 1290 1286 1291 1287 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)1288 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1293 1289 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) 1290 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1296 1291 1297 1292 DO jk = 1, jpkm1 - 1 … … 1316 1311 END DO 1317 1312 1318 ELSE IF (polynomial_type == 2) THEN ! Linear1313 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1319 1314 DO ji = 1, jpi 1320 1315 DO jj = 1, jpj … … 1347 1342 !! extrapolation is also permitted (no value limit) 1348 1343 !!---------------------------------------------------------------------- 1349 IMPLICIT NONE1350 1344 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1351 1345 REAL(wp) :: f ! result of the interpolation (extrapolation) 1352 1346 REAL(wp) :: zdeltx 1353 1347 !!---------------------------------------------------------------------- 1354 1348 ! 1355 1349 zdeltx = xr - xl 1356 IF( abs(zdeltx) <= 10._wp * EPSILON(x)) THEN1357 f = 0.5_wp * (fl + fr)1350 IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 1351 f = 0.5_wp * (fl + fr) 1358 1352 ELSE 1359 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx1353 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1360 1354 ENDIF 1361 1355 ! 1362 1356 END FUNCTION interp1 1363 1357 1364 1358 1365 FUNCTION interp2( x, a, b, c, d) RESULT(f)1359 FUNCTION interp2( x, a, b, c, d ) RESULT(f) 1366 1360 !!---------------------------------------------------------------------- 1367 1361 !! *** ROUTINE interp1 *** … … 1372 1366 !! 1373 1367 !!---------------------------------------------------------------------- 1374 IMPLICIT NONE1375 1368 REAL(wp), INTENT(in) :: x, a, b, c, d 1376 1369 REAL(wp) :: f ! value from the interpolation 1377 1370 !!---------------------------------------------------------------------- 1378 1371 ! 1379 1372 f = a + x* ( b + x * ( c + d * x ) ) 1380 1373 ! 1381 1374 END FUNCTION interp2 1382 1375 1383 1376 1384 FUNCTION interp3( x, a, b, c, d) RESULT(f)1377 FUNCTION interp3( x, a, b, c, d ) RESULT(f) 1385 1378 !!---------------------------------------------------------------------- 1386 1379 !! *** ROUTINE interp1 *** … … 1392 1385 !! 1393 1386 !!---------------------------------------------------------------------- 1394 IMPLICIT NONE1395 1387 REAL(wp), INTENT(in) :: x, a, b, c, d 1396 1388 REAL(wp) :: f ! value from the interpolation 1397 1389 !!---------------------------------------------------------------------- 1398 1390 ! 1399 1391 f = b + x * ( 2._wp * c + 3._wp * d * x) 1400 1392 ! 1401 1393 END FUNCTION interp3 1402 1394 1403 1395 1404 FUNCTION integ_spline( xl, xr, a, b, c, d) RESULT(f)1396 FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) 1405 1397 !!---------------------------------------------------------------------- 1406 1398 !! *** ROUTINE interp1 *** … … 1411 1403 !! 1412 1404 !!---------------------------------------------------------------------- 1413 IMPLICIT NONE1414 1405 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1415 1406 REAL(wp) :: za1, za2, za3 1416 1407 REAL(wp) :: f ! integration result 1417 1408 !!---------------------------------------------------------------------- 1418 1409 ! 1419 1410 za1 = 0.5_wp * b 1420 1411 za2 = c / 3.0_wp 1421 1412 za3 = 0.25_wp * d 1422 1413 ! 1423 1414 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1424 1415 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1425 1416 ! 1426 1417 END FUNCTION integ_spline 1427 1418
Note: See TracChangeset
for help on using the changeset viewer.
