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 3524 for branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 – NEMO

Ignore:
Timestamp:
2012-11-02T07:13:40+01:00 (11 years ago)
Author:
gm
Message:

Branch: dev_r3385_NOCS04_HAMF; #665. add USE lib_fortran when SIGN is used (TOP,OPA,LIM2&3) ; salt flux names start with sfx_ in LIM3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3523 r3524  
    1313   !!   'key_lim3'                                      LIM3 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_lat_acr    : lateral accretion of ice 
    16    !!---------------------------------------------------------------------- 
    17    USE par_oce          ! ocean parameters 
    18    USE dom_oce          ! domain variables 
    19    USE phycst           ! physical constants 
    20    USE sbc_oce          ! Surface boundary condition: ocean fields 
    21    USE sbc_ice          ! Surface boundary condition: ice fields 
    22    USE thd_ice          ! LIM thermodynamics 
    23    USE dom_ice          ! LIM domain 
    24    USE par_ice          ! LIM parameters 
    25    USE ice              ! LIM variables 
    26    USE limtab           ! LIM 2D <==> 1D 
    27    USE limcons          ! LIM conservation 
    28    USE in_out_manager   ! I/O manager 
    29    USE lib_mpp          ! MPP library 
    30    USE wrk_nemo         ! work arrays 
     15   !!   lim_lat_acr   : lateral accretion of ice 
     16   !!---------------------------------------------------------------------- 
     17   USE par_oce        ! ocean parameters 
     18   USE dom_oce        ! domain variables 
     19   USE phycst         ! physical constants 
     20   USE sbc_oce        ! Surface boundary condition: ocean fields 
     21   USE sbc_ice        ! Surface boundary condition: ice fields 
     22   USE thd_ice        ! LIM thermodynamics 
     23   USE dom_ice        ! LIM domain 
     24   USE par_ice        ! LIM parameters 
     25   USE ice            ! LIM variables 
     26   USE limtab         ! LIM 2D <==> 1D 
     27   USE limcons        ! LIM conservation 
     28   USE in_out_manager ! I/O manager 
     29   USE lib_mpp        ! MPP library 
     30   USE wrk_nemo       ! work arrays 
     31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3132 
    3233   IMPLICIT NONE 
     
    4546 
    4647   !!---------------------------------------------------------------------- 
    47    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     48   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4849   !! $Id$ 
    4950   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7778      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)       
    7879      !!------------------------------------------------------------------------ 
    79       INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices 
    80       INTEGER ::   layer, nbpac     ! local integers  
    81       INTEGER ::   zji, zjj, iter   !   -       - 
    82       REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde  ! local scalars 
     80      INTEGER  ::   ji,jj,jk,jl,jm   ! dummy loop indices 
     81      INTEGER  ::   layer, nbpac     ! local integers  
     82      INTEGER  ::   zji, zjj, iter   !   -       - 
     83      REAL(wp) ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde   ! local scalars 
    8384      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8485      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
     86      REAL(wp) ::   zcoef                                                       !   -      - 
    8587      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8688      CHARACTER (len = 15) :: fieldid 
     
    205207                  !------------- 
    206208                  ! C-grid wind stress components 
    207                   ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  ) & 
    208                      &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) / 2.0 
    209                   ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1) & 
    210                      &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) / 2.0 
     209                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   & 
     210                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp 
     211                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   & 
     212                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp 
    211213                  ! Square root of wind stress 
    212214                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
     
    222224                  !------------------- 
    223225                  ! C-grid ice velocity 
    224                   zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) )) 
    225                   zvgx  = zindb * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  ) & 
    226                      + u_ice(ji,jj    ) * tmu(ji  ,jj  ) ) / 2.0 
    227                   zvgy  = zindb * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1) & 
    228                      + v_ice(ji,jj    ) * tmv(ji  ,jj  ) ) / 2.0 
     226                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
     227                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    & 
     228                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp 
     229                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    & 
     230                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp 
    229231 
    230232                  !----------------------------------- 
     
    232234                  !----------------------------------- 
    233235                  ! absolute relative velocity 
    234                   zvrel2        = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + & 
    235                      ( zvfry - zvgy ) * ( zvfry - zvgy )   & 
    236                      , 0.15 * 0.15 ) 
    237                   zvrel(ji,jj)  = SQRT(zvrel2) 
     236                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
     237                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 
     238                  zvrel(ji,jj)  = SQRT( zvrel2 ) 
    238239 
    239240                  !--------------------- 
     
    241242                  !--------------------- 
    242243                  hicol(ji,jj) = zhicrit + 0.1  
    243                   hicol(ji,jj) = zhicrit + hicol(ji,jj) /      &  
    244                      ( hicol(ji,jj) * hicol(ji,jj) - & 
    245                      zhicrit * zhicrit ) * ztwogp * zvrel2 
     244                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    & 
     245                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
     246 
     247!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1) 
     248!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01 
     249!!gm                therefore the 2 lines with hicol can be replaced by 1 line: 
     250!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2 
     251!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop 
    246252 
    247253                  iter = 1 
     
    278284      DO jj = 1, jpj 
    279285         DO ji = 1, jpi 
    280             IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN 
     286            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN 
    281287               nbpac = nbpac + 1 
    282288               npac( nbpac ) = (jj - 1) * jpi + ji 
    283                IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 
    284                   jiindex_1d = nbpac 
    285                ENDIF 
     289               IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac 
    286290            ENDIF 
    287291         END DO 
    288292      END DO 
    289293 
    290       IF( ln_nicep ) THEN 
    291          WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    292       ENDIF 
     294      IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    293295 
    294296      !------------------------------ 
     
    302304         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) ) 
    303305         DO jl = 1, jpl 
    304             CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) ) 
    305             CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) ) 
    306             CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , jpi, jpj, npac(1:nbpac) ) 
     306            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     307            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     308            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    307309            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    308310            DO jk = 1, nlay_i 
     
    311313         END DO ! jl 
    312314 
    313          CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif , jpi, jpj, npac(1:nbpac) ) 
    314          CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif , jpi, jpj, npac(1:nbpac) ) 
    315          CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  , jpi, jpj, npac(1:nbpac) ) 
    316          CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv , jpi, jpj, npac(1:nbpac) ) 
     315         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) ) 
     316         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) ) 
     317         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) ) 
     318         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) ) 
    317319         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) ) 
    318          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol , jpi, jpj, npac(1:nbpac) ) 
    319          CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel , jpi, jpj, npac(1:nbpac) ) 
     320         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
     321         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
    320322 
    321323         !------------------------------------------------------------------------------! 
     
    392394         ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 
    393395         DO ji = 1, nbpac 
    394             fseqv_1d  (ji) = fseqv_1d  (ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
     396            sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
    395397            rdm_ice_1d(ji) = rdm_ice_1d(ji) +                 rhoic * zv_newice(ji) 
    396398         END DO ! ji 
     
    400402         !------------------------------------ 
    401403         DO ji = 1, nbpac 
    402             ! Volume 
    403             zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    404             zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    405             vt_i_init(zji,zjj)   = vt_i_init(zji,zjj) + zv_newice(ji) 
    406             ! Energy 
    407             zde                  = ze_newice(ji) / unit_fac 
    408             zde                  = zde * area(zji,zjj) * zv_newice(ji) 
    409             et_i_init(zji,zjj)   = et_i_init(zji,zjj) + zde 
     404            zji = MOD( npac(ji) - 1 , jpi ) + 1 
     405            zjj =    ( npac(ji) - 1 ) / jpi + 1 
     406            ! 
     407            zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji) 
     408            ! 
     409            vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)             ! volume 
     410            et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde                       ! Energy 
     411 
    410412         END DO 
    411413 
     
    417419         !----------------- 
    418420         DO ji = 1, nbpac 
    419             za_newice(ji)     = zv_newice(ji) / zh_newice(ji) 
    420             ! diagnostic 
    421             zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    422             zjj                  = ( npac(ji) - 1 ) / jpi + 1 
     421            zji = MOD( npac(ji) - 1 , jpi ) + 1 
     422            zjj =    ( npac(ji) - 1 ) / jpi + 1 
     423            za_newice(ji) = zv_newice(ji) / zh_newice(ji) 
    423424            diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 
    424425         END DO !ji 
     
    438439         !------------------------------------------- 
    439440         ! If lateral ice growth gives an ice concentration gt 1, then 
    440          ! we keep the excessive volume in memory and attribute it later 
    441          ! to bottom accretion 
    442          DO ji = 1, nbpac 
    443             ! vectorize 
    444             IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN 
    445                zda_res(ji)    = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 
    446                zdv_res(ji)    = zda_res(ji) * zh_newice(ji)  
    447                za_newice(ji)  = za_newice(ji) - zda_res(ji) 
    448                zv_newice(ji)  = zv_newice(ji) - zdv_res(ji) 
     441         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
     442         DO ji = 1, nbpac 
     443            IF ( za_newice(ji)  >  ( 1._wp - zat_i_ac(ji) ) ) THEN 
     444               zda_res(ji)   = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 
     445               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
     446               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     447               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji) 
    449448            ELSE 
    450                zda_res(ji) = 0.0 
    451                zdv_res(ji) = 0.0 
     449               zda_res(ji) = 0._wp 
     450               zdv_res(ji) = 0._wp 
    452451            ENDIF 
    453452         END DO ! ji 
     
    459458         DO jl = 1, jpl 
    460459            DO ji = 1, nbpac 
    461                IF(  hi_max   (jl-1)  <  zh_newice(ji)   .AND.   & 
    462                   & zh_newice(ji)    <= hi_max   (jl)         ) THEN 
     460               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   & 
     461                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN 
    463462                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 
    464463                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 
     
    466465                  zcatac  (ji)    = jl 
    467466               ENDIF 
    468             END DO ! ji 
    469          END DO ! jl 
     467            END DO 
     468         END DO 
    470469 
    471470         !---------------------------------- 
     
    483482            DO ji = 1, nbpac 
    484483               jl = zcatac(ji) 
    485                zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 
     484               zqold   = ze_i_ac(ji,jk,jl)      ! [ J.m-3 ] 
    486485               zalphai = MIN( zhice_old(ji,jl) *   jk       / nlay_i , zh_newice(ji) )   & 
    487486                  &    - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) ) 
     
    489488                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / nlay_i   & 
    490489                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       & 
    491                   + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i ) 
     490                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i ) 
    492491            END DO 
    493492         END DO 
     
    529528               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    & 
    530529                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce 
    531                zdummy(ji,jl)    = zv_i_ac(ji,jl)/MAX(za_i_ac(ji,jl),epsi10)*zindb      ! thickness of residual ice 
     530               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice 
    532531            END DO 
    533532         END DO 
     
    613612            END DO 
    614613         END DO 
    615          CALL tab_1d_2d( nbpac, fseqv  , npac(1:nbpac), fseqv_1d  (1:nbpac), jpi, jpj ) 
     614         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 
    616615         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 
    617616         ! 
     
    623622      DO jl = 1, jpl 
    624623         DO jk = 1, nlay_i          ! heat content in 10^9 Joules 
    625             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i  / unit_fac  
     624            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac  
    626625         END DO 
    627626      END DO 
Note: See TracChangeset for help on using the changeset viewer.