Changeset 5905 for branches/2015/dev_r5151_UKMO_ISF
- Timestamp:
- 2015-11-20T17:59:57+01:00 (8 years ago)
- Location:
- branches/2015/dev_r5151_UKMO_ISF/NEMOGCM
- Files:
-
- 2 edited
-
CONFIG/SHARED/namelist_ref (modified) (2 diffs)
-
NEMO/OPA_SRC/SBC/sbcisf.F90 (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/CONFIG/SHARED/namelist_ref
r5834 r5905 440 440 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! 441 441 ! nn_isf == 4 442 sn_qisf = 'rnfisf' , -12 ,'sohflisf', .false. , .true. , 'yearly' , '' , '' 443 sn_fwfisf = 'rnfisf' , -12 ,'sowflisf', .false. , .true. , 'yearly' , '' , '' 442 sn_fwfisf = 'rnfisf' , -12 ,'sowflisf', .false. , .true. , 'yearly' , '' , '' 444 443 ! nn_isf == 3 445 sn_rnfisf = 'runoffs' , -12 ,'sofwfisf', .false., .true. , 'yearly' , '' , ''444 sn_rnfisf = 'rnfisf' , -12 ,'sofwfisf', .false. , .true. , 'yearly' , '' , '' 446 445 ! nn_isf == 2 and 3 447 sn_depmax_isf = 'r unoffs' , -12 ,'sozisfmax' , .false., .true. , 'yearly' , '' , ''448 sn_depmin_isf = 'r unoffs' , -12 ,'sozisfmin' , .false., .true. , 'yearly' , '' , ''446 sn_depmax_isf = 'rnfisf' , -12 ,'sozisfmax' , .false. , .true. , 'yearly' , '' , '' 447 sn_depmin_isf = 'rnfisf' , -12 ,'sozisfmin' , .false. , .true. , 'yearly' , '' , '' 449 448 ! nn_isf == 2 450 sn_Leff_isf = 'rnfisf' , 0 ,'Leff' , .false., .true. , 'yearly' , '' , ''449 sn_Leff_isf = 'rnfisf' , -12 ,'Leff' , .false. , .true. , 'yearly' , '' , '' 451 450 ! for all case 452 451 nn_isf = 1 ! ice shelf melting/freezing … … 458 457 rn_gammas0 = 1.0e-4 ! gammas coefficient used in blk formula 459 458 ! only for nn_isf = 1 or 4 460 rn_hisf_tbl = 30. ! thickness of the top boundary layer 459 rn_hisf_tbl = 30. ! thickness of the top boundary layer (Losh et al. 2008) 461 460 ! 0 => thickness of the tbl = thickness of the first wet cell 462 461 ! only for nn_isf = 1 463 nn_isfblk = 1 ! 1 ISOMIP ; 2 : 3 equation formulation, Jenkins et al. 1991 462 nn_isfblk = 1 ! 1 ISOMIP like: 2 equations formulation (Hunter et al., 2006) 463 ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 464 464 nn_gammablk = 1 ! 0 = cst Gammat (= gammat/s) 465 465 ! 1 = velocity dependend Gamma (u* * gammat/s) (Jenkins et al. 2010) 466 ! if you want to keep the cd as in global config, adjust rn_gammat0 to compensate 467 ! 2 = velocity and stability dependent Gamma Holland et al. 1999 466 ! 2 = velocity and stability dependent Gamma (Holland et al. 1999) 468 467 / 469 468 !----------------------------------------------------------------------- -
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5834 r5905 68 68 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 69 69 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 70 TYPE(FLD_N) , PUBLIC :: sn_ qisf, sn_fwfisf !: information about the runoff file to be read71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qisf,sf_fwfisf70 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf file to be read 71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 72 72 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 73 73 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf … … 85 85 86 86 SUBROUTINE sbc_isf(kt) 87 INTEGER, INTENT(in) :: kt ! ocean time step 88 INTEGER :: ji, jj, jk, inum, ierror 89 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 87 INTEGER, INTENT(in) :: kt ! ocean time step 88 INTEGER :: ji, jj, jk ! loop index 89 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 90 INTEGER :: inum, ierror 91 INTEGER :: ios ! Local integer output status for namelist read 90 92 REAL(wp) :: zhk 91 REAL(wp) :: zt_frz, zpress93 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 92 94 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 93 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 94 INTEGER :: ios ! Local integer output status for namelist read 95 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 95 96 ! 96 97 !!--------------------------------------------------------------------- 97 98 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf , & 98 & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 99 ! 99 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 100 ! 101 ! allocation 102 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 100 103 ! 101 104 ! ! ====================== ! … … 148 151 CALL iom_close(inum) 149 152 ! 150 risfLeff = risfLeff*1000 !: convertion in m153 risfLeff = risfLeff*1000.0_wp !: convertion in m 151 154 END IF 152 155 … … 179 182 180 183 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 181 ALLOCATE( sf_fwfisf(1), sf_qisf(1),STAT=ierror )184 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 182 185 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 183 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )184 186 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 185 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' )186 187 END IF 187 188 … … 197 198 198 199 ! determine the deepest level influenced by the boundary layer 199 ! test on tmask useless ????? 200 DO jk = ikt, mbkt(ji,jj) 200 DO jk = ikt+1, mbkt(ji,jj) 201 201 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 202 202 END DO … … 247 247 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 248 248 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 249 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! f resh waterflux from the isf (fwfisf <0 mean melting)250 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat 249 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 250 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 251 251 stbl(:,:) = soce 252 252 … … 254 254 ! specified fwf and heat flux forcing beneath the ice shelf 255 255 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 256 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 257 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 258 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 259 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux 256 fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 257 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 260 258 stbl(:,:) = soce 261 262 259 END IF 260 263 261 ! compute tsc due to isf 264 ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable). 265 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 266 zt_frz = -1.9_wp !CALL eos_fzp( tsn(ji,jj,jk,jp_sal), zt_frz, zpress ) 267 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz * r1_rau0 ! 262 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 263 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 264 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 265 DO jj = 1,jpj 266 DO ji = 1,jpi 267 zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj)) 268 END DO 269 END DO 270 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 268 271 269 ! salt effect already take into account in vertical advection272 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 270 273 risf_tsc(:,:,jp_sal) = 0.0_wp 271 274 … … 276 279 CALL lbc_lnk(qisf(:,:) ,'T',1.) 277 280 278 IF( kt == nit000 ) THEN 281 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 279 282 IF( ln_rstart .AND. & ! Restart: read in restart file 280 283 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN … … 293 296 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 294 297 END IF 298 299 ! deallocation 300 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 295 301 296 302 END SUBROUTINE sbc_isf … … 360 366 ! after verif with UNESCO, wrong sign in BG eq. 2 361 367 ! Calculate freezing temperature 362 zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 363 CALL eos_fzp(tsb(ji,jj,jk,jp_sal), zt_frz, zpress) 368 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 364 369 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 365 370 ENDDO … … 398 403 INTEGER, INTENT(in) :: kt ! ocean time step 399 404 ! 400 REAL(wp), DIMENSION(:,:), POINTER :: zfrz ,zpress,zti405 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 401 406 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 402 407 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b … … 411 416 !!--------------------------------------------------------------------- 412 417 ! 413 ! coeficient for linearisation of tfreez 414 zlamb1=-0.0575_wp 415 zlamb2= 0.0901_wp 416 zlamb3=-7.61e-04 418 ! coeficient for linearisation of potential tfreez 419 ! Crude approximation for pressure (but commonly used) 420 zlamb1=-0.0573_wp 421 zlamb2= 0.0832_wp 422 zlamb3=-7.53e-08 * grav * rau0 417 423 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 418 424 ! 419 CALL wrk_alloc( jpi,jpj, zfrz , z press, zti , zgammat, zgammas)420 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b)425 CALL wrk_alloc( jpi,jpj, zfrz , zgammat, zgammas ) 426 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 421 427 422 428 ! initialisation 423 zpress (:,:)=0.0_wp424 429 zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0; 425 430 zhtflx (:,:)=0.0_wp ; zhtflx_b(:,:)=0.0_wp ; 426 431 zfwflx (:,:)=0.0_wp 427 ! 428 ! Crude approximation for pressure (but commonly used) 429 ! 1e-04 to convert from Pa to dBar 430 DO jj = 1, jpj 431 DO ji = 1, jpi 432 zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 433 END DO 434 END DO 435 436 ! Calculate in-situ temperature (ref to surface) 437 zti(:,:) = tinsitu( ttbl, stbl, zpress ) 438 432 433 ! compute ice shelf melting 439 434 nit = 1 ; lit = .TRUE. 440 435 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 441 IF (nn_isfblk == 1) THEN ! ISOMIP case 436 IF (nn_isfblk == 1) THEN 437 ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 442 438 ! Calculate freezing temperature 443 CALL eos_fzp( s ss_m(:,:), zfrz(:,:), zpress(:,:) )439 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 444 440 445 441 ! compute gammat every where (2d) … … 449 445 DO jj = 1, jpj 450 446 DO ji = 1, jpi 451 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*( zti(ji,jj)-zfrz(ji,jj))447 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 452 448 zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf 453 449 END DO … … 455 451 456 452 ! Compute heat flux and upward fresh water flux 457 ! For genuine ISOMIP protocol this should probably be something like 458 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 459 fwfisf(:,:) = zfwflx(:,:) * ( soce / MAX(stbl(:,:),zeps)) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 453 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 454 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 460 455 461 456 ELSE IF (nn_isfblk == 2 ) THEN 462 ! More complicated 3 equation thermodynamics as in MITgcm457 ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 463 458 ! compute gammat every where (2d) 464 459 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 465 460 466 461 ! compute upward heat flux zhtflx and upward water flux zwflx 462 ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 467 463 DO jj = 1, jpj 468 464 DO ji = 1, jpi … … 472 468 zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 473 469 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 474 zeps6=zeps4- zti(ji,jj)470 zeps6=zeps4-ttbl(ji,jj) 475 471 zeps7=zeps4-tsurf 476 472 zaqe=zlamb1 * (zeps1 + zeps3) … … 485 481 IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 486 482 487 ! compute t freeze 483 ! compute t freeze (eq. 22) 488 484 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 489 485 490 486 ! zfwflx is upward water flux 491 487 ! zhtflx is upward heat flux (out of ocean) 492 ! compute the upward water and heat flux 488 ! compute the upward water and heat flux (eq. 28 and eq. 29) 493 489 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 494 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * ( zti(ji,jj) - zfrz(ji,jj) )490 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 495 491 END DO 496 492 END DO … … 523 519 IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 524 520 ! 525 CALL wrk_dealloc( jpi,jpj, zfrz , z press, zti , zgammat, zgammas)526 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b)521 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas ) 522 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 527 523 ! 528 524 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') … … 563 559 IF ( nn_gammablk == 0 ) THEN 564 560 !! gamma is constant (specified in namelist) 561 !! ISOMIP formulation (Hunter et al, 2006) 565 562 gt(:,:) = rn_gammat0 566 563 gs(:,:) = rn_gammas0 … … 569 566 !! gamma is assume to be proportional to u* 570 567 !! Jenkins et al., 2010, JPO, p2298-2312 571 572 !! compute ustar 568 !! Adopted by Asay-Davis et al. (2015) 569 570 !! compute ustar (eq. 24) 573 571 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 574 572 … … 604 602 zcoef = 0.5 / fse3w(ji,jj,ikt) 605 603 ! ! shear of horizontal velocity 606 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) 604 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 607 605 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 608 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) 606 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 609 607 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 610 608 ! ! richardson number (minimum value set to zero) … … 664 662 665 663 REAL(wp) :: ze3, zhk 666 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl 667 668 INTEGER :: ji,jj,jk 669 INTEGER :: ikt,ikb 670 664 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 665 666 INTEGER :: ji,jj,jk ! loop index 667 INTEGER :: ikt,ikb ! top and bottom index of the tbl 668 669 ! allocation 671 670 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 672 673 ! get first and last level of tbl671 672 ! initialisation 674 673 varout(:,:)=0._wp 674 675 ! compute U in the top boundary layer at T- point 675 676 IF (cptin == 'U') THEN 676 677 DO jj = 1,jpj … … 681 682 682 683 ! determine the deepest level influenced by the boundary layer 683 ! test on tmask useless ????? 684 DO jk = ikt, mbku(ji,jj) 684 DO jk = ikt+1, mbku(ji,jj) 685 685 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 686 686 END DO … … 706 706 END IF 707 707 708 ! compute V in the top boundary layer at T- point 708 709 IF (cptin == 'V') THEN 709 710 DO jj = 1,jpj … … 714 715 715 716 ! determine the deepest level influenced by the boundary layer 716 ! test on tmask useless ????? 717 DO jk = ikt, mbkv(ji,jj) 717 DO jk = ikt+1, mbkv(ji,jj) 718 718 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 719 719 END DO … … 739 739 END IF 740 740 741 ! compute T in the top boundary layer at T- point 741 742 IF (cptin == 'T') THEN 742 743 DO jj = 1,jpj … … 757 758 END DO 758 759 END IF 760 761 ! mask mean tbl value 759 762 varout(:,:) = varout(:,:) * ssmask(:,:) 760 763 764 ! deallocation 761 765 CALL wrk_dealloc( jpi,jpj, zhisf_tbl ) 762 766 … … 795 799 796 800 ! determine the deepest level influenced by the boundary layer 797 DO jk = ikt , mbkt(ji,jj)801 DO jk = ikt+1, mbkt(ji,jj) 798 802 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 799 803 END DO … … 825 829 ! 826 830 END SUBROUTINE sbc_isf_div 827 828 FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 829 !!---------------------------------------------------------------------- 830 !! *** ROUTINE eos_init *** 831 !! 832 !! ** Purpose : Compute the in-situ temperature [Celcius] 833 !! 834 !! ** Method : 835 !! 836 !! Reference : Bryden,h.,1973,deep-sea res.,20,401-408 837 !!---------------------------------------------------------------------- 838 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptem ! potential temperature [Celcius] 839 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 840 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ppress ! pressure [dBar] 841 REAL(wp), DIMENSION(:,:), POINTER :: pti ! in-situ temperature [Celcius] 842 REAL(wp) :: zt, zs, zp, zh, zq, zxk 843 INTEGER :: ji, jj ! dummy loop indices 844 ! 845 CALL wrk_alloc( jpi,jpj, pti ) 846 ! 847 DO jj=1,jpj 848 DO ji=1,jpi 849 zh = ppress(ji,jj) 850 ! Theta1 851 zt = ptem(ji,jj) 852 zs = psal(ji,jj) 853 zp = 0.0 854 zxk= zh * fsatg( zs, zt, zp ) 855 zt = zt + 0.5 * zxk 856 zq = zxk 857 ! Theta2 858 zp = zp + 0.5 * zh 859 zxk= zh*fsatg( zs, zt, zp ) 860 zt = zt + 0.29289322 * ( zxk - zq ) 861 zq = 0.58578644 * zxk + 0.121320344 * zq 862 ! Theta3 863 zxk= zh * fsatg( zs, zt, zp ) 864 zt = zt + 1.707106781 * ( zxk - zq ) 865 zq = 3.414213562 * zxk - 4.121320344 * zq 866 ! Theta4 867 zp = zp + 0.5 * zh 868 zxk= zh * fsatg( zs, zt, zp ) 869 pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0 870 END DO 871 END DO 872 ! 873 CALL wrk_dealloc( jpi,jpj, pti ) 874 ! 875 END FUNCTION tinsitu 876 ! 877 FUNCTION fsatg( pfps, pfpt, pfphp ) 878 !!---------------------------------------------------------------------- 879 !! *** FUNCTION fsatg *** 880 !! 881 !! ** Purpose : Compute the Adiabatic laspse rate [Celcius].[decibar]^-1 882 !! 883 !! ** Reference : Bryden,h.,1973,deep-sea res.,20,401-408 884 !! 885 !! ** units : pressure pfphp decibars 886 !! temperature pfpt deg celsius (ipts-68) 887 !! salinity pfps (ipss-78) 888 !! adiabatic fsatg deg. c/decibar 889 !!---------------------------------------------------------------------- 890 REAL(wp) :: pfps, pfpt, pfphp 891 REAL(wp) :: fsatg 892 ! 893 fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp & 894 & +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt & 895 & +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp & 896 & +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.) & 897 & +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 898 ! 899 END FUNCTION fsatg 900 !!====================================================================== 831 !!====================================================================== 901 832 END MODULE sbcisf
Note: See TracChangeset
for help on using the changeset viewer.
