New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4946 for branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 – NEMO

Ignore:
Timestamp:
2014-12-02T10:38:20+01:00 (9 years ago)
Author:
cetlod
Message:

2014/dev_MERGE_2014 : merge in changes from dev_CNRS_CICE

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r4938 r4946  
    5454   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2 
    5555   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 
    56    INTEGER(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
     56#ifdef key_agrif 
     57   ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals 
     58   REAL(wp),    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
    5759                                                                                          !: (first wet level and last level include in the tbl) 
     60#else 
     61   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
     62#endif 
     63 
    5864 
    5965   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
     
    303309      sbc_isf_alloc = 0       ! set to zero if no array to be allocated 
    304310      IF( .NOT. ALLOCATED( qisf ) ) THEN 
    305          ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj), fwfisf(jpi,jpj), & 
    306                &    fwfisf_b(jpi,jpj), misfkt(jpi,jpj), rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj),     & 
    307                &    rzisf_tbl(jpi,jpj), misfkb(jpi,jpj), ttbl(jpi,jpj), stbl(jpi,jpj), utbl(jpi,jpj), & 
    308                &    vtbl(jpi, jpj), risfLeff(jpi,jpj), rhisf_tbl_0(jpi,jpj), ralpha(jpi,jpj), STAT= sbc_isf_alloc ) 
     311         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts)              , & 
     312               &    qisf(jpi,jpj)     , fwfisf(jpi,jpj)     , fwfisf_b(jpi,jpj)   , & 
     313               &    rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
     314               &    ttbl(jpi,jpj)     , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
     315               &    vtbl(jpi, jpj)    , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
     316               &    ralpha(jpi,jpj)   , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
     317               &    STAT= sbc_isf_alloc ) 
    309318         ! 
    310319         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc ) 
     
    363372             ! Calculate freezing temperature 
    364373                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04  
    365                 zt_frz = tfreez1D(tsb(ji,jj,ik,jp_sal), zpress)  
     374                zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress)  
    366375                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp 
    367376             ENDDO 
     
    445454      zti(:,:)=tinsitu( ttbl, stbl, zpress ) 
    446455! Calculate freezing temperature 
    447       zfrz(:,:)=tfreez( sss_m(:,:), zpress ) 
     456      zfrz(:,:)=eos_fzp( sss_m(:,:), zpress ) 
    448457 
    449458       
     
    526535                     nit = nit + 1 
    527536                     IF (nit .GE. 51) THEN 
    528                         WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 
     537                        WRITE(numout,*) "sbcisf : too many iteration ... ", & 
     538                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 
    529539                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    530540                     END IF 
     
    584594      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence  
    585595      REAL(wp) :: zcoef                      ! temporary coef 
    586       REAL(wp) :: zrhos, zalbet, zbeta, zthermal, zhalin 
    587       REAL(wp) :: zt, zs, zh 
     596      REAL(wp) :: zdep 
    588597      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant 
    589598      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 
    590599      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 
    591600      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s 
     601      REAL(wp), DIMENSION(2) :: zts, zab 
    592602      !!--------------------------------------------------------------------- 
    593603      ! 
     
    656666 
    657667      !! compute bouyancy  
    658                IF( nn_eos < 1) THEN 
    659                   zt     = ttbl(ji,jj) 
    660                   zs     = stbl(ji,jj) - 35.0 
    661                   zh     = fsdepw(ji,jj,ikt) 
    662                   !  potential volumic mass 
    663                   zrhos  = rhop(ji,jj,ikt) 
    664                   zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta 
    665                      &                               - 0.203814e-03 ) * zt   & 
    666                      &                               + 0.170907e-01 ) * zt   & 
    667                      &   + 0.665157e-01                                      & 
    668                      &   +     ( - 0.678662e-05 * zs                         & 
    669                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    670                      &   +   ( ( - 0.302285e-13 * zh                         & 
    671                      &           - 0.251520e-11 * zs                         & 
    672                      &           + 0.512857e-12 * zt * zt           ) * zh   & 
    673                      &           - 0.164759e-06 * zs                         & 
    674                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    675                      &                               + 0.380374e-04 ) * zh 
    676  
    677                   zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta 
    678                      &                            - 0.301985e-05 ) * zt      & 
    679                      &   + 0.785567e-03                                      & 
    680                      &   + (     0.515032e-08 * zs                           & 
    681                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    682                      &   +(  (   0.121551e-17 * zh                           & 
    683                      &         - 0.602281e-15 * zs                           & 
    684                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    685                      &                             + 0.408195e-10   * zs     & 
    686                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    687                      &                             - 0.121555e-07 ) * zh 
    688  
    689                   zthermal = zbeta * zalbet / ( rcp * zrhos + epsln ) 
    690                   zhalin   = zbeta * stbl(ji,jj) * rcs 
    691                ELSE 
    692                   zrhos    = rhop(ji,jj,ikt) + rau0 * ( 1. - tmask(ji,jj,ikt) ) 
    693                   zthermal = rn_alpha / ( rcp * zrhos + epsln ) 
    694                   zhalin   = rn_beta * stbl(ji,jj) * rcs 
    695                ENDIF 
     668               zts(jp_tem) = ttbl(ji,jj) 
     669               zts(jp_sal) = stbl(ji,jj) 
     670               zdep        = fsdepw(ji,jj,ikt) 
     671               ! 
     672               CALL eos_rab( zts, zdep, zab ) 
     673                  ! 
    696674      !! compute length scale  
    697                zbuofdep = grav * ( zthermal * zqhisf - zhalin * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     675               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    698676 
    699677      !! compute Monin Obukov Length 
     
    766744               ! level partially include in ice shelf boundary layer  
    767745               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    768                IF (cptin == 'T') varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
    769                IF (cptin == 'U') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 
    770                IF (cptin == 'V') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 
     746               IF (cptin == 'T') & 
     747                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     748               IF (cptin == 'U') & 
     749                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 
     750               IF (cptin == 'V') & 
     751                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 
    771752            END IF 
    772753         END DO 
     
    796777      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    797778      !! 
    798       INTEGER(wp)  ::   ji, jj, jk   ! dummy loop indices 
    799       INTEGER(wp)  ::   ikt, ikb  
    800       INTEGER(wp)  ::   nk_isf 
     779      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     780      INTEGER  ::   ikt, ikb  
     781      INTEGER  ::   nk_isf 
    801782      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl 
    802783      REAL(wp)     ::   zfact     ! local scalar 
     
    834815               ! level fully include in the ice shelf boundary layer 
    835816               DO jk = ikt, ikb - 1 
    836                   phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
     817                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
     818                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    837819               END DO 
    838820               ! level partially include in ice shelf boundary layer  
    839                phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
     821               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 
     822                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    840823            !==   ice shelf melting mass distributed over several levels   ==! 
    841824         END DO 
Note: See TracChangeset for help on using the changeset viewer.