Ignore:
Timestamp:
2015-11-20T17:59:57+01:00 (5 years ago)
Author:
mathiot
Message:

ISF: based on N. Jourdain comments, compute fwf from potential temp instead of insitu, extra comment + rm useless line

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r5834 r5905  
    6868!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    6969   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 read 
    71    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf 
     70   TYPE(FLD_N)       , PUBLIC ::   sn_fwfisf            !: information about the isf file to be read 
     71   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_fwfisf 
    7272   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read 
    7373   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf            
     
    8585  
    8686  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 
    9092    REAL(wp)            :: zhk 
    91     REAL(wp)            ::   zt_frz, zpress 
     93    REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
    9294    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 
    9596      ! 
    9697      !!--------------------------------------------------------------------- 
    9798      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  ) 
    100103      ! 
    101104      !                                         ! ====================== ! 
     
    148151               CALL iom_close(inum) 
    149152               ! 
    150                risfLeff = risfLeff*1000           !: convertion in m 
     153               risfLeff = risfLeff*1000.0_wp           !: convertion in m 
    151154            END IF 
    152155 
     
    179182             
    180183            ! 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 ) 
    182185            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) ) 
    184186            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' ) 
    186187         END IF 
    187188          
     
    197198 
    198199               ! 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) 
    201201                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    202202               END DO 
     
    247247            ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    248248            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    249             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
    250             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
     249            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf flux from the isf (fwfisf <0 mean melting)  
     250            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat flux 
    251251            stbl(:,:)   = soce 
    252252 
     
    254254            ! specified fwf and heat flux forcing beneath the ice shelf 
    255255            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 
    260258            stbl(:,:)   = soce 
    261  
    262259         END IF 
     260 
    263261         ! 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(:,:) ) 
    268271          
    269          ! salt effect already take into account in vertical advection 
     272         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 
    270273         risf_tsc(:,:,jp_sal) = 0.0_wp 
    271274 
     
    276279         CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
    277280 
    278          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     281         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    279282            IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
    280283                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 
     
    293296         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf) 
    294297      END IF 
     298 
     299      ! deallocation 
     300      CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 
    295301   
    296302  END SUBROUTINE sbc_isf 
     
    360366             ! after verif with UNESCO, wrong sign in BG eq. 2 
    361367             ! 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)  
    364369                zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
    365370             ENDDO 
     
    398403      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    399404      ! 
    400       REAL(wp), DIMENSION(:,:), POINTER ::   zfrz,zpress,zti 
     405      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz 
    401406      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
    402407      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
     
    411416      !!--------------------------------------------------------------------- 
    412417      ! 
    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 
    417423      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav') 
    418424      ! 
    419       CALL wrk_alloc( jpi,jpj, zfrz  , zpress, 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 ) 
    421427 
    422428      ! initialisation 
    423       zpress (:,:)=0.0_wp 
    424429      zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0; 
    425430      zhtflx (:,:)=0.0_wp     ; zhtflx_b(:,:)=0.0_wp    ; 
    426431      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 
    439434      nit = 1 ; lit = .TRUE. 
    440435      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) 
    442438            ! Calculate freezing temperature 
    443             CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress(:,:) ) 
     439            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 
    444440 
    445441            ! compute gammat every where (2d) 
     
    449445            DO jj = 1, jpj 
    450446               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)) 
    452448                  zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf 
    453449               END DO 
     
    455451 
    456452            ! 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(:,:) 
    460455 
    461456         ELSE IF (nn_isfblk == 2 ) THEN 
    462             ! More complicated 3 equation thermodynamics as in MITgcm 
     457            ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)  
    463458            ! compute gammat every where (2d) 
    464459            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
    465460 
    466461            ! 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) 
    467463            DO jj = 1, jpj 
    468464               DO ji = 1, jpi 
     
    472468                  zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 
    473469                  zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    474                   zeps6=zeps4-zti(ji,jj) 
     470                  zeps6=zeps4-ttbl(ji,jj) 
    475471                  zeps7=zeps4-tsurf 
    476472                  zaqe=zlamb1 * (zeps1 + zeps3) 
     
    485481                  IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    486482 
    487                   ! compute t freeze 
     483                  ! compute t freeze (eq. 22) 
    488484                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    489485   
    490486                  ! zfwflx is upward water flux 
    491487                  ! 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) 
    493489                  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) )  
    495491               END DO 
    496492            END DO 
     
    523519      IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 
    524520      !  
    525       CALL wrk_dealloc( jpi,jpj, zfrz  , zpress, 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 ) 
    527523      ! 
    528524      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav') 
     
    563559      IF      ( nn_gammablk == 0 ) THEN 
    564560      !! gamma is constant (specified in namelist) 
     561      !! ISOMIP formulation (Hunter et al, 2006) 
    565562         gt(:,:) = rn_gammat0 
    566563         gs(:,:) = rn_gammas0 
     
    569566      !! gamma is assume to be proportional to u*  
    570567      !! 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) 
    573571         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
    574572 
     
    604602                  zcoef = 0.5 / fse3w(ji,jj,ikt) 
    605603                  !                                            ! 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  )  & 
    607605                     &             -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  )  & 
    609607                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
    610608                  !                                            ! richardson number (minimum value set to zero) 
     
    664662 
    665663      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 
    671670      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    672  
    673       ! get first and last level of tbl 
     671       
     672      ! initialisation 
    674673      varout(:,:)=0._wp 
     674    
     675      ! compute U in the top boundary layer at T- point  
    675676      IF (cptin == 'U') THEN 
    676677         DO jj = 1,jpj 
     
    681682 
    682683            ! 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) 
    685685                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
    686686               END DO 
     
    706706      END IF 
    707707 
     708      ! compute V in the top boundary layer at T- point  
    708709      IF (cptin == 'V') THEN 
    709710         DO jj = 1,jpj 
     
    714715 
    715716            ! 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) 
    718718                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
    719719               END DO 
     
    739739      END IF 
    740740 
     741      ! compute T in the top boundary layer at T- point  
    741742      IF (cptin == 'T') THEN 
    742743         DO jj = 1,jpj 
     
    757758         END DO 
    758759      END IF 
     760 
     761      ! mask mean tbl value 
    759762      varout(:,:) = varout(:,:) * ssmask(:,:) 
    760763 
     764      ! deallocation 
    761765      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )       
    762766 
     
    795799 
    796800               ! determine the deepest level influenced by the boundary layer 
    797                DO jk = ikt, mbkt(ji,jj) 
     801               DO jk = ikt+1, mbkt(ji,jj) 
    798802                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    799803               END DO 
     
    825829      ! 
    826830   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   !!====================================================================== 
    901832END MODULE sbcisf 
Note: See TracChangeset for help on using the changeset viewer.