- Timestamp:
- 2015-04-09T20:32:14+02:00 (9 years ago)
- Location:
- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 edited
-
SBC/sbcisf.F90 (modified) (16 diffs)
-
TRA/zpshde.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5120 r5204 87 87 88 88 SUBROUTINE sbc_isf(kt) 89 INTEGER, INTENT(in) :: kt ! ocean time step 90 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 91 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 92 REAL(wp) :: rmin 93 REAL(wp) :: zhk 94 CHARACTER(len=256) :: cfisf, cvarzisf, cvarhisf ! name for isf file 95 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 96 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 97 INTEGER :: ios ! Local integer output status for namelist read 89 INTEGER, INTENT(in) :: kt ! ocean time step 90 INTEGER :: ji, jj, jk, inum, ierror 91 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 92 REAL(wp) :: zhk 93 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 94 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 95 INTEGER :: ios ! Local integer output status for namelist read 98 96 ! 99 97 !!--------------------------------------------------------------------- … … 125 123 IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf 126 124 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 125 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 126 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 127 127 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 128 128 IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM) … … 150 150 !: read effective lenght (BG03) 151 151 IF (nn_isf == 2) THEN 152 ! Read Data and save some integral values152 cvarLeff = 'soLeff' 153 153 CALL iom_open( sn_Leff_isf%clname, inum ) 154 cvarLeff = 'soLeff' !: variable name for Efficient Length scale155 154 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 156 155 CALL iom_close(inum) … … 237 236 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 238 237 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 239 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')240 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')238 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 239 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 241 240 ! iom print 242 241 CALL iom_put('ttbl',ttbl(:,:)) … … 296 295 ! 297 296 ! output 298 CALL iom_put('qisf' , qisf)299 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl (:,:)/ soce )297 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 298 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl / soce ) 300 299 END IF 301 300 … … 343 342 INTEGER, INTENT ( in ) :: kt 344 343 345 INTEGER :: ji, jj, jk, jish !temporary integer 346 INTEGER :: ijkmin 347 INTEGER :: ii, ij, ik 348 INTEGER :: inum 349 350 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 351 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 352 REAL(wp) :: zt_frz ! freezing point temperature at depth z 353 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 344 INTEGER :: ji, jj, jk !temporary integer 345 346 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 347 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 348 REAL(wp) :: zt_frz ! freezing point temperature at depth z 349 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 354 350 355 351 !!---------------------------------------------------------------------- … … 360 356 DO ji = 1, jpi 361 357 DO jj = 1, jpj 362 ik = misfkt(ji,jj)358 jk = misfkt(ji,jj) 363 359 !! Initialize arrays to 0 (each step) 364 360 zt_sum = 0.e0_wp 365 IF ( ik .GT. 1 ) THEN366 ! 3. -----------the average temperature between 200m and 600m ---------------------361 IF ( jk .GT. 1 ) THEN 362 ! 1. -----------the average temperature between 200m and 600m --------------------- 367 363 DO jk = misfkt(ji,jj),misfkb(ji,jj) 368 364 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 369 365 ! after verif with UNESCO, wrong sign in BG eq. 2 370 366 ! Calculate freezing temperature 371 zpress = grav*rau0*fsdept(ji,jj, ik)*1.e-04372 zt_frz = eos_fzp(tsb(ji,jj, ik,jp_sal), zpress)373 zt_sum = zt_sum + (tsn(ji,jj, ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik) ! sum temp367 zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 368 zt_frz = eos_fzp(tsb(ji,jj,jk,jp_sal), zpress) 369 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 374 370 ENDDO 375 371 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 376 372 377 ! 4. ------------Net heat flux and fresh water flux due to the ice shelf373 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 378 374 ! For those corresponding to zonal boundary 379 375 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 380 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj, ik)376 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 381 377 382 378 fwfisf(ji,jj) = qisf(ji,jj) / lfusisf !fresh water flux kg/(m2s) … … 407 403 INTEGER, INTENT(in) :: kt ! ocean time step 408 404 ! 409 LOGICAL :: ln_isomip = .true. 410 REAL(wp), DIMENSION(:,:), POINTER :: zfrz,zpress,zti 411 REAL(wp), DIMENSION(:,:), POINTER :: zgammat2d, zgammas2d 412 !REAL(wp), DIMENSION(:,:), POINTER :: zqisf, zfwfisf 405 REAL(wp), DIMENSION(:,:), POINTER :: zfrz,zpress,zti 406 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 407 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 413 408 REAL(wp) :: zlamb1, zlamb2, zlamb3 414 409 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 415 410 REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 416 REAL(wp) :: zfwflx, zhtflx, zhtflx_b 417 REAL(wp) :: zgammat, zgammas 418 REAL(wp) :: zeps = -1.e-20_wp !== Local constant initialization ==! 411 REAL(wp) :: zeps = 1.e-20_wp 412 REAL(wp) :: zerr 419 413 INTEGER :: ji, jj ! dummy loop indices 420 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers421 INTEGER :: ierror ! return error code422 LOGICAL :: lit=.TRUE.423 414 INTEGER :: nit 415 LOGICAL :: lit 424 416 !!--------------------------------------------------------------------- 425 417 ! 426 418 ! coeficient for linearisation of tfreez 427 zlamb1=-0.0575 428 zlamb2= 0.0901419 zlamb1=-0.0575_wp 420 zlamb2= 0.0901_wp 429 421 zlamb3=-7.61e-04 430 422 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 431 423 ! 432 CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 433 434 zcfac=0.0_wp 435 IF (ln_conserve) zcfac=1.0_wp 436 zpress(:,:)=0.0_wp 437 zgammat2d(:,:)=0.0_wp 438 zgammas2d(:,:)=0.0_wp 439 ! 440 ! 424 CALL wrk_alloc( jpi,jpj, zfrz , zpress, zti , zgammat, zgammas ) 425 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b ) 426 427 ! initialisation 428 IF (ln_conserve) THEN ; zcfac=1.0_wp 429 ELSE ; zcfac=0.0_wp 430 ENDIF 431 zpress (:,:)=0.0_wp 432 zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0; 433 zhtflx (:,:)=0.0_wp ; zhtflx_b(:,:)=0.0_wp ; 434 zfwflx (:,:)=0.0_wp 435 ! 436 ! Crude approximation for pressure (but commonly used) 437 ! 1e-04 to convert from Pa to dBar 441 438 !CDIR COLLAPSE 442 439 DO jj = 1, jpj 443 440 DO ji = 1, jpi 444 ! Crude approximation for pressure (but commonly used)445 ! 1e-04 to convert from Pa to dBar446 441 zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 447 !448 442 END DO 449 443 END DO 450 444 451 ! Calculate in-situ temperature (ref to surface) 452 zti(:,:)=tinsitu( ttbl, stbl, zpress ) 453 ! Calculate freezing temperature 454 zfrz(:,:)=eos_fzp( sss_m(:,:), zpress ) 455 456 457 zhtflx=0._wp ; zfwflx=0._wp 458 IF (nn_isfblk == 1) THEN 459 DO jj = 1, jpj 460 DO ji = 1, jpi 461 IF (mikt(ji,jj) > 1 ) THEN 462 nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp 463 DO WHILE ( lit ) 464 ! compute gamma 465 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 466 ! zhtflx is upward heat flux (out of ocean) 467 zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 468 ! zwflx is upward water flux 469 zfwflx = - zhtflx/lfusisf 470 ! test convergence and compute gammat 471 IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 472 473 nit = nit + 1 474 IF (nit .GE. 100) THEN 475 !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj 476 !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx 477 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 478 END IF 479 ! save gammat and compute zhtflx_b 480 zgammat2d(ji,jj)=zgammat 481 zhtflx_b = zhtflx 482 END DO 483 484 qisf(ji,jj) = - zhtflx 485 ! For genuine ISOMIP protocol this should probably be something like 486 fwfisf(ji,jj) = zfwflx * ( soce / MAX(stbl(ji,jj),zeps)) 487 ELSE 488 fwfisf(ji,jj) = 0._wp 489 qisf(ji,jj) = 0._wp 490 END IF 491 ! 445 ! Calculate in-situ temperature (ref to surface) 446 zti(:,:) = tinsitu( ttbl, stbl, zpress ) 447 448 ! Calculate freezing temperature 449 zfrz(:,:) = eos_fzp( sss_m, zpress ) 450 451 nit = 1 ; lit = .TRUE. 452 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 453 IF (nn_isfblk == 1) THEN ! ISOMIP case 454 ! compute gammat every where (2d) 455 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 456 457 ! compute upward heat flux zhtflx and upward water flux zwflx 458 DO jj = 1, jpj 459 DO ji = 1, jpi 460 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 461 zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf 462 END DO 492 463 END DO 493 END DO 494 495 ELSE IF (nn_isfblk == 2 ) THEN 496 497 ! More complicated 3 equation thermodynamics as in MITgcm 498 !CDIR COLLAPSE 499 DO jj = 2, jpj 500 DO ji = 2, jpi 501 IF (mikt(ji,jj) > 1 ) THEN 502 nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp 503 DO WHILE ( lit ) 504 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 505 506 zeps1=rcp*rau0*zgammat 507 zeps2=lfusisf*rau0*zgammas 508 zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj) 509 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 510 zeps6=zeps4-zti(ji,jj) 511 zeps7=zeps4-tsurf 512 zaqe=zlamb1 * (zeps1 + zeps3) 513 zaqer=0.5/zaqe 514 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 515 zcqe=zeps2*stbl(ji,jj) 516 zdis=zbqe*zbqe-4.0*zaqe*zcqe 517 ! Presumably zdis can never be negative because gammas is very small compared to gammat 518 zsfrz=(-zbqe-SQRT(zdis))*zaqer 519 IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer 520 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 464 465 ! Compute heat flux and upward fresh water flux 466 ! For genuine ISOMIP protocol this should probably be something like 467 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 468 fwfisf(:,:) = zfwflx(:,:) * ( soce / MAX(stbl(:,:),zeps)) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 469 470 ELSE IF (nn_isfblk == 2 ) THEN 471 ! More complicated 3 equation thermodynamics as in MITgcm 472 ! compute gammat every where (2d) 473 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 474 475 ! compute upward heat flux zhtflx and upward water flux zwflx 476 DO jj = 1, jpj 477 DO ji = 1, jpi 478 ! compute coeficient to solve the 2nd order equation 479 zeps1=rcp*rau0*zgammat(ji,jj) 480 zeps2=lfusisf*rau0*zgammas(ji,jj) 481 zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 482 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 483 zeps6=zeps4-zti(ji,jj) 484 zeps7=zeps4-tsurf 485 zaqe=zlamb1 * (zeps1 + zeps3) 486 zaqer=0.5/MIN(zaqe,-zeps) 487 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 488 zcqe=zeps2*stbl(ji,jj) 489 zdis=zbqe*zbqe-4.0*zaqe*zcqe 490 491 ! Presumably zdis can never be negative because gammas is very small compared to gammat 492 ! compute s freeze 493 IF (zsfrz .GE. 0.0_wp) THEN ; zsfrz=(-zbqe-SQRT(zdis))*zaqer 494 ELSE ; zsfrz=(-zbqe+SQRT(zdis))*zaqer 495 ENDIF 496 497 ! compute t freeze 498 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 521 499 522 ! zfwflx is upward water flux 523 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz ) 524 ! zhtflx is upward heat flux (out of ocean) 500 ! zfwflx is upward water flux 501 ! zhtflx is upward heat flux (out of ocean) 525 502 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 526 zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) )527 ! zwflx is upward water flux528 503 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 529 zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 530 ! test convergence and compute gammat 531 IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 532 533 nit = nit + 1 534 IF (nit .GE. 51) THEN 535 WRITE(numout,*) "sbcisf : too many iteration ... ", & 536 & zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 537 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 538 END IF 539 ! save gammat and compute zhtflx_b 540 zgammat2d(ji,jj)=zgammat 541 zgammas2d(ji,jj)=zgammas 542 zhtflx_b = zhtflx 543 544 END DO 545 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 546 qisf(ji,jj) = - zhtflx 547 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 548 fwfisf(ji,jj) = zfwflx 549 ELSE 550 fwfisf(ji,jj) = 0._wp 551 qisf(ji,jj) = 0._wp 552 ENDIF 553 ! 504 ! compute the upward water and heat flux 505 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 506 zhtflx(ji,jj) = ( zgammat(ji,jj) * rau0 - zcfac * zfwflx(ji,jj) ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 507 zfwflx(ji,jj) = ( zgammas(ji,jj) * rau0 - zcfac * zfwflx(ji,jj) ) * (zsfrz - stbl(ji,jj)) / MAX(stbl(ji,jj),zeps) 508 END DO 554 509 END DO 555 END DO 556 ENDIF 557 ! lbclnk 558 CALL lbc_lnk(zgammas2d(:,:),'T',1.) 559 CALL lbc_lnk(zgammat2d(:,:),'T',1.) 510 511 ! compute heat and water flux 512 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 513 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 514 515 ENDIF 516 517 ! define if we need to iterate 518 IF ( nn_gammablk .LT. 2 ) THEN ; lit = .FALSE. 519 ELSE 520 ! check total number of iteration 521 IF (nit .GE. 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 522 ELSE ; nit = nit + 1 523 ENDIF 524 525 ! compute error between 2 iterations 526 ! if needed save gammat and compute zhtflx_b for next iteration 527 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 528 IF ( zerr .LE. 0.01 ) THEN ; lit = .FALSE. 529 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 530 ENDIF 531 END IF 532 END DO 533 ! 560 534 ! output 561 CALL iom_put('isfgammat', zgammat2d)562 CALL iom_put('isfgammas', zgammas2d)563 !564 !CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zqisf, zfwfisf)565 CALL wrk_dealloc( jpi,jpj, zf rz,zpress,zti, zgammat2d, zgammas2d)535 IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 536 IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 537 ! 538 CALL wrk_dealloc( jpi,jpj, zfrz , zpress, zti , zgammat, zgammas ) 539 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx, zhtflx_b ) 566 540 ! 567 541 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') … … 569 543 END SUBROUTINE sbc_isf_cav 570 544 571 SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf , ji, jj, lit)545 SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf ) 572 546 !!---------------------------------------------------------------------- 573 547 !! ** Purpose : compute the coefficient echange for heat flux … … 578 552 !! Jenkins et al., 2010, JPO, p2298-2312 579 553 !!--------------------------------------------------------------------- 580 REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf 581 INTEGER , INTENT(in) :: ji,jj 582 LOGICAL , INTENT(inout) :: lit 554 REAL(wp), DIMENSION(:,:), INTENT(out) :: gt, gs 555 REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf 583 556 584 557 INTEGER :: ikt ! loop index 585 REAL(wp) :: zut, zvt, zustar ! U, V at T point and friction velocity 558 INTEGER :: ji,jj 559 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity 586 560 REAL(wp) :: zdku, zdkv ! U, V shear 587 561 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number … … 596 570 REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 597 571 REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 598 REAL(wp) :: rcs = 1.0e-3_wp ! conversion: mm/s ==> m/s572 REAL(wp) :: zeps = 1.0e-20_wp ! conversion: mm/s ==> m/s 599 573 REAL(wp), DIMENSION(2) :: zts, zab 600 574 !!--------------------------------------------------------------------- 601 ! 602 IF( nn_gammablk == 0 ) THEN 575 CALL wrk_alloc( jpi,jpj, zustar ) 576 ! 577 IF ( nn_gammablk == 0 ) THEN 603 578 !! gamma is constant (specified in namelist) 604 gt = rn_gammat0605 gs = rn_gammas0606 lit = .FALSE. 579 gt(:,:) = rn_gammat0 580 gs(:,:) = rn_gammas0 581 607 582 ELSE IF ( nn_gammablk == 1 ) THEN 608 583 !! gamma is assume to be proportional to u* 609 !! WARNING in case of Losh 2008 tbl parametrization,610 !! you have to used the mean value of u in the boundary layer)611 !! not yet coded612 584 !! Jenkins et al., 2010, JPO, p2298-2312 613 ikt = mikt(ji,jj)614 !! Compute U and V at T points615 ! zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) )616 ! zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) )617 zut = utbl(ji,jj)618 zvt = vtbl(ji,jj)619 585 620 586 !! compute ustar 621 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 622 !! Compute mean value over the TBL 587 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) ) 623 588 624 589 !! Compute gammats 625 gt = zustar* rn_gammat0626 gs = zustar* rn_gammas0627 lit = .FALSE.590 gt(:,:) = zustar(:,:) * rn_gammat0 591 gs(:,:) = zustar(:,:) * rn_gammas0 592 628 593 ELSE IF ( nn_gammablk == 2 ) THEN 629 594 !! gamma depends of stability of boundary layer 630 !! WARNING in case of Losh 2008 tbl parametrization,631 !! you have to used the mean value of u in the boundary layer)632 !! not yet coded633 595 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 634 596 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 597 !! compute ustar 598 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) ) 599 600 !! compute Pr and Sc number (can be improved) 601 zPr = 13.8 602 zSc = 2432.0 603 604 !! compute gamma mole 605 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 606 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 607 608 !! compute gamma 609 DO ji=2,jpi 610 DO jj=2,jpj 635 611 ikt = mikt(ji,jj) 636 612 637 !! Compute U and V at T points 638 zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) ) 639 zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) ) 640 641 !! compute ustar 642 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 643 IF (zustar == 0._wp) THEN ! only for kt = 1 I think 644 gt = rn_gammat0 645 gs = rn_gammas0 613 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think 614 gt = rn_gammat0 615 gs = rn_gammas0 646 616 ELSE 647 617 !! compute Rc number (as done in zdfric.F90) 648 zcoef = 0.5 / fse3w(ji,jj,ikt) 649 ! ! shear of horizontal velocity 650 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 651 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 652 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 653 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 654 ! ! richardson number (minimum value set to zero) 655 zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 656 657 !! compute Pr and Sc number (can be improved) 658 zPr = 13.8 659 zSc = 2432.0 660 661 !! compute gamma mole 662 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 663 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 618 zcoef = 0.5 / fse3w(ji,jj,ikt) 619 ! ! shear of horizontal velocity 620 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 621 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 622 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 623 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 624 ! ! richardson number (minimum value set to zero) 625 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 664 626 665 627 !! compute bouyancy 666 zts(jp_tem) = ttbl(ji,jj)667 zts(jp_sal) = stbl(ji,jj)668 zdep = fsdepw(ji,jj,ikt)669 !670 CALL eos_rab( zts, zdep, zab )628 zts(jp_tem) = ttbl(ji,jj) 629 zts(jp_sal) = stbl(ji,jj) 630 zdep = fsdepw(ji,jj,ikt) 631 ! 632 CALL eos_rab( zts, zdep, zab ) 671 633 ! 672 634 !! compute length scale 673 zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf) !!!!!!!!!!!!!!!!!!!!!!!!!!!!635 zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 674 636 675 637 !! compute Monin Obukov Length 676 ! Maximum boundary layer depth677 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001678 ! Compute Monin obukhov length scale at the surface and Ekman depth:679 zmob = zustar** 3 / (vkarmn * (zbuofdep + epsln))680 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)638 ! Maximum boundary layer depth 639 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 640 ! Compute Monin obukhov length scale at the surface and Ekman depth: 641 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + epsln)) 642 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 681 643 682 644 !! compute eta* (stability parameter) 683 zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))645 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 684 646 685 647 !! compute the sublayer thickness 686 zhnu = 5 * znu / zustar 648 zhnu = 5 * znu / zustar(ji,jj) 649 687 650 !! compute gamma turb 688 zgturb = 1/vkarmn * LOG(zustar* zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &689 & + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn651 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 652 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 690 653 691 654 !! compute gammats 692 gt = zustar/ (zgturb + zgmolet)693 gs = zustar/ (zgturb + zgmoles)655 gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 656 gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 694 657 END IF 658 END DO 659 END DO 660 CALL lbc_lnk(gt(:,:),'T',1.) 661 CALL lbc_lnk(gs(:,:),'T',1.) 695 662 END IF 663 CALL wrk_dealloc( jpi,jpj, zustar ) 696 664 697 665 END SUBROUTINE … … 701 669 !! *** SUBROUTINE sbc_isf_tbl *** 702 670 !! 703 !! ** Purpose : compute mean T/S/U/V in the boundary layer 671 !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 704 672 !! 705 673 !!---------------------------------------------------------------------- … … 726 694 DO jj = 2,jpj 727 695 DO ji = 2,jpi 728 IF (ssmask(ji,jj) == 1) THEN 729 ikt = mkt(ji,jj) 730 ikb = mkb(ji,jj) 731 732 ! level fully include in the ice shelf boundary layer 733 DO jk = ikt, ikb - 1 734 ze3 = fse3t_n(ji,jj,jk) 735 IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 736 IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 737 & * r1_hisf_tbl(ji,jj) * ze3 738 IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 739 & * r1_hisf_tbl(ji,jj) * ze3 740 END DO 741 742 ! level partially include in ice shelf boundary layer 743 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 744 IF (cptin == 'T') & 745 & varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 746 IF (cptin == 'U') & 747 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 748 IF (cptin == 'V') & 749 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 750 END IF 696 ikt = mkt(ji,jj) 697 ikb = mkb(ji,jj) 698 699 ! level fully include in the ice shelf boundary layer 700 DO jk = ikt, ikb - 1 701 ze3 = fse3t_n(ji,jj,jk) 702 IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 703 IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 704 & * r1_hisf_tbl(ji,jj) * ze3 705 IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 706 & * r1_hisf_tbl(ji,jj) * ze3 707 END DO 708 709 ! level partially include in ice shelf boundary layer 710 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 711 IF (cptin == 'T') & 712 & varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 713 IF (cptin == 'U') & 714 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 715 IF (cptin == 'V') & 716 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 751 717 END DO 752 718 END DO 719 varout(:,:) = varout(:,:) * ssmask(:,:) 753 720 754 721 CALL wrk_dealloc( jpi,jpj, mkt, mkb ) … … 777 744 INTEGER :: ji, jj, jk ! dummy loop indices 778 745 INTEGER :: ikt, ikb 779 INTEGER :: nk_isf 780 REAL(wp) :: zhk, z1_hisf_tbl, zhisf_tbl 781 REAL(wp) :: zfact ! local scalar 746 REAL(wp) :: zhk 747 REAL(wp) :: zfact ! local scalar 782 748 !!---------------------------------------------------------------------- 783 749 ! … … 795 761 ! test on tmask useless ????? 796 762 DO jk = ikt, mbkt(ji,jj) 797 !IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk763 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 798 764 END DO 799 765 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 839 805 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ppress ! pressure [dBar] 840 806 REAL(wp), DIMENSION(:,:), POINTER :: pti ! in-situ temperature [Celcius] 841 ! REAL(wp) :: fsatg842 ! REAL(wp) :: pfps, pfpt, pfphp843 807 REAL(wp) :: zt, zs, zp, zh, zq, zxk 844 808 INTEGER :: ji, jj ! dummy loop indices -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r5200 r5204 268 268 DO jj = 1, jpjm1 269 269 DO ji = 1, jpim1 270 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 271 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 270 271 iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 272 ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 272 273 ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 273 274 ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) … … 279 280 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 280 281 ! gradient of tracers 281 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )282 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 282 283 ELSE ! case 2 283 284 zmaxu = -ze3wu / fse3w(ji,jj,iku) … … 285 286 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 286 287 ! gradient of tracers 287 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )288 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 288 289 ENDIF 289 290 ! … … 294 295 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 295 296 ! gradient of tracers 296 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )297 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 297 298 ELSE ! case 2 298 299 zmaxv = -ze3wv / fse3w(ji,jj,ikv) … … 300 301 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 301 302 ! gradient of tracers 302 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 303 ENDIF 303 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 304 ENDIF 305 304 306 END DO 305 307 END DO … … 314 316 DO jj = 1, jpjm1 315 317 DO ji = 1, jpim1 318 316 319 iku = mbku(ji,jj) 317 320 ikv = mbkv(ji,jj) … … 337 340 DO jj = 1, jpjm1 338 341 DO ji = 1, jpim1 342 339 343 iku = mbku(ji,jj) 340 344 ikv = mbkv(ji,jj) … … 342 346 ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 343 347 344 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 345 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 346 ENDIF 347 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 348 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 349 ENDIF 350 END DO 351 END DO 348 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 349 ELSE ; pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 350 ENDIF 351 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 352 ELSE ; pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 353 ENDIF 354 355 END DO 356 END DO 357 352 358 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 353 359 ! … … 357 363 DO jj = 1, jpjm1 358 364 DO ji = 1, jpim1 359 iku = miku(ji,jj) ;ikup1 = miku(ji,jj) + 1360 ikv = mikv(ji,jj) ;ikvp1 = mikv(ji,jj) + 1365 iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 366 ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 361 367 ! 362 368 ! (ISF) case partial step top and bottom in adjacent cell in vertical … … 366 372 ze3wu = gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 367 373 ze3wv = gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 374 368 375 ! i- direction 369 376 IF( ze3wu >= 0._wp ) THEN ! case 1 370 zmaxu = ze3wu / fse3w(ji+1,jj,iku +1)371 ! interpolated values of tracers 372 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku +1,jn) - pta(ji+1,jj,iku,jn) )373 ! gradient of tracers 374 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )375 ELSE ! case 2 376 zmaxu = - ze3wu / fse3w(ji,jj,iku +1)377 ! interpolated values of tracers 378 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku +1,jn) - pta(ji,jj,iku,jn) )377 zmaxu = ze3wu / fse3w(ji+1,jj,ikup1) 378 ! interpolated values of tracers 379 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 380 ! gradient of tracers 381 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 382 ELSE ! case 2 383 zmaxu = - ze3wu / fse3w(ji,jj,ikup1) 384 ! interpolated values of tracers 385 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 379 386 ! gradient of tracers 380 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )387 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 381 388 ENDIF 382 389 ! 383 390 ! j- direction 384 391 IF( ze3wv >= 0._wp ) THEN ! case 1 385 zmaxv = ze3wv / fse3w(ji,jj+1,ikv+1) 386 ! interpolated values of tracers 387 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 388 ! gradient of tracers 389 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 390 ELSE ! case 2 391 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) 392 ! interpolated values of tracers 393 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 394 ! gradient of tracers 395 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 396 ENDIF 397 END DO!! 398 END DO!! 399 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 392 zmaxv = ze3wv / fse3w(ji,jj+1,ikvp1) 393 ! interpolated values of tracers 394 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 395 ! gradient of tracers 396 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 397 ELSE ! case 2 398 zmaxv = - ze3wv / fse3w(ji,jj,ikvp1) 399 ! interpolated values of tracers 400 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 401 ! gradient of tracers 402 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 403 ENDIF 404 405 END DO 406 END DO 407 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 400 408 ! 401 409 END DO … … 403 411 ! horizontal derivative of density anomalies (rd) 404 412 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 405 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 406 DO jj = 1, jpjm1 407 DO ji = 1, jpim1 413 pgrui(:,:) =0.0_wp; pgrvi(:,:) =0.0_wp; 414 DO jj = 1, jpjm1 415 DO ji = 1, jpim1 416 408 417 iku = miku(ji,jj) 409 418 ikv = mikv(ji,jj) … … 430 439 DO jj = 1, jpjm1 431 440 DO ji = 1, jpim1 432 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 433 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 441 442 iku = miku(ji,jj) 443 ikv = mikv(ji,jj) 434 444 ze3wu = gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 435 445 ze3wv = gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv) 436 446 437 IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj ) - prd(ji,jj,iku) )! i: 1438 ELSE ; pgrui(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj ,iku) - zri(ji,jj ) )! i: 2439 ENDIF 440 441 IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji ,jj ) - prd(ji,jj,ikv) ) ! j: 1442 ELSE ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji ,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2443 ENDIF 444 445 END DO 446 END DO 447 CALL lbc_lnk( pgrui , 'U', -1. ) ;CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions447 IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 448 ELSE ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj ,iku) - zri(ji,jj ) ) ! i: 2 449 ENDIF 450 451 IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji ,jj ) - prd(ji,jj,ikv) ) ! j: 1 452 ELSE ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji ,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 453 ENDIF 454 455 END DO 456 END DO 457 CALL lbc_lnk( pgrui , 'U', -1. ); CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 448 458 ! 449 459 END IF
Note: See TracChangeset
for help on using the changeset viewer.
