Changeset 5078


Ignore:
Timestamp:
2015-02-11T16:15:11+01:00 (6 years ago)
Author:
clem
Message:

LIM3: cosmetic changes to increase readability and performance

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5076 r5078  
    165165 
    166166   !                                     !!** ice-thickness distribution namelist (namiceitd) ** 
    167    INTEGER , PUBLIC ::   jpl              !: number of ice  categories  
    168    INTEGER , PUBLIC ::   nlay_i           !: number of ice  layers  
    169    INTEGER , PUBLIC ::   nlay_s           !: number of snow layers  
    170167   INTEGER , PUBLIC ::   nn_catbnd        !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 
    171168   REAL(wp), PUBLIC ::   rn_himean        !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 
     
    223220   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( rn_ecc * rn_ecc ) 
    224221   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cio 
    225  
     222   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i 
     223   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s  
     224   ! 
    226225   !                                     !!** switch for presence of ice or not  
    227226   REAL(wp), PUBLIC ::   rswitch 
    228  
     227   ! 
    229228   !                                     !!** define some parameters  
    230229   REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number  
     
    368367   !!-------------------------------------------------------------------------- 
    369368   !                                                  !!: ** Namelist namicerun read in sbc_lim_init ** 
    370    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
    371    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
    372    LOGICAL               , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F) 
    373    LOGICAL               , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F) 
    374    REAL(wp)              , PUBLIC ::   rn_amax         !: maximum ice concentration 
     369   INTEGER          , PUBLIC ::   jpl             !: number of ice  categories  
     370   INTEGER          , PUBLIC ::   nlay_i          !: number of ice  layers  
     371   INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers  
     372   CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
     373   CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
     374   LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F) 
     375   LOGICAL          , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F) 
     376   REAL(wp)         , PUBLIC ::   rn_amax         !: maximum ice concentration 
    375377   ! 
    376378   !!-------------------------------------------------------------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5076 r5078  
    210210            &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw  
    211211 
    212          zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
     212         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic ) 
    213213 
    214214         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5067 r5078  
    326326               ! recompute ht_i, ht_s avoiding out of bounds values 
    327327               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
    328                ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
     328               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
    329329 
    330330               ! ice volume, salt content, age content 
     
    347347 
    348348                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    349                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s 
     349                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    350350               END DO ! ji 
    351351            END DO ! jj 
     
    368368 
    369369                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    370                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / nlay_i 
     370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
    371371               END DO ! ji 
    372372            END DO ! jj 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5070 r5078  
    11241124               ! clem: if sst>0, then ersw <0 (is that possible?) 
    11251125               zsstK  = sst_m(ji,jj) + rt0 
    1126                ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 
     1126               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i 
    11271127 
    11281128               ! heat flux to the ocean 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5076 r5078  
    337337      DO jl = 1, jpl 
    338338         DO jk = 1, nlay_i 
    339             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i ) 
     339            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i 
    340340         END DO 
    341341      END DO 
     
    347347      DO jl = 1, jpl 
    348348         DO jk = 1, nlay_s 
    349             e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s ) 
     349            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s 
    350350         END DO 
    351351      END DO 
     
    479479            ! Conversion q(S,T) -> T (second order equation) 
    480480            zaaa          =  cpic 
    481             zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus 
     481            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) * r1_rhoic - lfus 
    482482            zccc          =  lfus * ( ztmelts - rtt ) 
    483483            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5070 r5078  
    145145      DO jk = 1, nlay_i 
    146146         DO ji = kideb, kiut 
    147             h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     147            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    148148            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 
    149149         ENDDO 
     
    188188      ! 
    189189      DO ji = kideb, kiut      
    190          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     190         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    191191      END DO 
    192192      ! 
     
    199199      DO jk = 1, nlay_i 
    200200         DO ji = kideb, kiut 
    201             zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 
     201            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 
    202202            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 
    203203         END DO 
     
    228228         ! thickness change 
    229229         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji)  
    230          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 
     230         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 
    231231         ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    232232         zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     
    255255         zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )       
    256256         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 
    257          zh_s  (ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     257         zh_s  (ji) = ht_s_1d(ji) * r1_nlay_s 
    258258 
    259259         ENDIF 
     
    311311      DO ji = kideb, kiut 
    312312         dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    313          zh_s(ji)       = ht_s_1d(ji) / REAL( nlay_s ) 
     313         zh_s(ji)       = ht_s_1d(ji) * r1_nlay_s 
    314314      END DO ! ji 
    315315 
     
    335335      DO jk = 1, nlay_i 
    336336         DO ji = kideb, kiut  
    337             zEi            = - q_i_1d(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0] 
     337            zEi            = - q_i_1d(ji,jk) * r1_rhoic             ! Specific enthalpy of layer k [J/kg, <0] 
    338338 
    339339            ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K] 
     
    345345            zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    346346 
    347             zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0] 
     347            zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
    348348 
    349349            zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     
    506506               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
    507507 
    508                   zEi               = - q_i_1d(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0) 
     508                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    509509 
    510510                  !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
     
    535535               ELSE                               !!! Basal melting 
    536536 
    537                   zEi               = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    538  
    539                   zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 
    540  
    541                   zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0) 
    542  
    543                   zfmdt             = - zq_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0) 
    544  
    545                   zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change 
     537                  zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     538 
     539                  zEw               = rcp * ( ztmelts - rtt )    ! Specific enthalpy of meltwater (J/kg, <0) 
     540 
     541                  zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     542 
     543                  zfmdt             = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     544 
     545                  zdeltah(ji,jk)    = - zfmdt * r1_rhoic         ! Gross thickness change 
    546546 
    547547                  zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     
    627627         ! Salinity of snow ice 
    628628         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    629          zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
     629         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
    630630 
    631631         ! entrapment during snow ice formation 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5076 r5078  
    179179      zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
    180180      DO ji = kideb, kiut 
    181          zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
    182             &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) )  
     181         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     182            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )  
    183183      END DO 
    184184 
     
    189189         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not 
    190190         ! layer thickness 
    191          zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 
    192          zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 
     191         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i 
     192         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    193193      END DO 
    194194 
     
    202202      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    203203         DO ji = kideb , kiut 
    204             z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 
     204            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 
    205205         END DO 
    206206      END DO 
     
    208208      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    209209         DO ji = kideb , kiut 
    210             z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 
     210            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 
    211211         END DO 
    212212      END DO 
     
    576576               !------------------------------------------------------------------------------| 
    577577               ! 
    578                IF (t_su_1d(ji) < rtt) THEN 
     578               IF ( t_su_1d(ji) < rtt ) THEN 
    579579                  ! 
    580580                  !------------------------------------------------------------------------------| 
     
    600600                  !!case of only one layer in the ice (surface & ice equations are altered) 
    601601 
    602                   IF (nlay_i.eq.1) THEN 
     602                  IF ( nlay_i == 1 ) THEN 
    603603                     ztrid(ji,numeqmin(ji),1)    =  0.0 
    604604                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0 
     
    631631 
    632632                  !!case of only one layer in the ice (surface & ice equations are altered) 
    633                   IF (nlay_i.eq.1) THEN 
     633                  IF ( nlay_i == 1 ) THEN 
    634634                     ztrid(ji,numeqmin(ji),1)  =  0.0 
    635635                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     
    760760      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    761761      DO ji = kideb, kiut 
    762          zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  & 
    763             &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
     762         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & 
     763            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
    764764         IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC 
    765765            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice  
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r5064 r5078  
    102102      ! new layer thickesses 
    103103      DO ji = kideb, kiut 
    104          zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i )   
     104         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i   
    105105      ENDDO 
    106106 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5070 r5078  
    362362         DO ji = 1, nbpac 
    363363 
    364             zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg] 
    365  
    366             zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg] 
     364            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg] 
     365 
     366            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg] 
    367367                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)  
    368368                                                                    
     
    371371            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)  
    372372                                                                   ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point    
    373             zv_newice(ji) = - zfmdt / rhoic 
     373            zv_newice(ji) = - zfmdt * r1_rhoic 
    374374 
    375375            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux   
     
    458458            DO jk = 1, nlay_i 
    459459               DO ji = 1, nbpac 
    460                   h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 
     460                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 
    461461                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 
    462462               END DO 
     
    525525               DO ji = 1, jpi 
    526526                  ! heat content in J/m2 
    527                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp )  
     527                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i  
    528528               END DO 
    529529            END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5067 r5078  
    196196                  ! 
    197197                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation) 
    198                   zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus 
     198                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i * r1_rhoic - lfus 
    199199                  zccc       =  lfus * (ztmelts-rtt) 
    200200                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 
     
    280280      !!------------------------------------------------------------------ 
    281281      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index 
    282       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal 
    283       REAL(wp) ::   zswi0, zswi01, zswibal, zargtemp , zs_zero    
     282      REAL(wp) ::   zfac0, zfac1, zsal 
     283      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero    
    284284      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha 
    285285      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
     
    311311         END DO 
    312312         ! 
    313          dummy_fac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf 
    314          dummy_fac1 = zsi1 / ( zsi1 - zsi0 ) 
     313         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf 
     314         zfac1 = zsi1 / ( zsi1 - zsi0 ) 
    315315         ! 
    316316         zalpha(:,:,:) = 0._wp 
     
    322322                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws  
    323323                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) )  
    324                   ! If 2.sm_i GE sss_m then zswibal = 1 
     324                  ! If 2.sm_i GE sss_m then rswitch = 1 
    325325                  ! this is to force a constant salinity profile in the Baltic Sea 
    326                   zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
    327                   zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
    328                   zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal ) 
    329                END DO 
    330             END DO 
    331          END DO 
    332  
    333          dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile 
     326                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
     327                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 
     328                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 
     329               END DO 
     330            END DO 
     331         END DO 
     332 
     333         ! Computation of the profile 
    334334         DO jl = 1, jpl 
    335335            DO jk = 1, nlay_i 
     
    337337                  DO ji = 1, jpi 
    338338                     !                                      ! linear profile with 0 at the surface 
    339                      zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac 
     339                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 
    340340                     !                                      ! weighting the profile 
    341341                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 
     
    357357         DO jl = 1, jpl 
    358358            DO jk = 1, nlay_i 
    359                zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
     359               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 
    360360               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  ) 
    361361               s_i(:,:,jk,jl) =  zsal 
     
    387387                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    388388                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    389                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 
     389                     &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) 
    390390               END DO 
    391391            END DO 
     
    417417                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  ) 
    418418                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   & 
    419                      &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp) 
     419                     &                   * v_i(ji,jj,jl) * r1_nlay_i 
    420420                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    421421                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 ) 
     
    439439      INTEGER  ::   ji, jk    ! dummy loop indices 
    440440      INTEGER  ::   ii, ij    ! local integers 
    441       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars 
    442       REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      - 
     441      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars 
     442      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      - 
    443443      ! 
    444444      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s 
     
    466466         ! Weighting factor between zs_zero and zs_inf 
    467467         !--------------------------------------------- 
    468          dummy_fac0 = 1._wp / ( zsi0 - zsi1 ) 
    469          dummy_fac1 = zsi1 / ( zsi1 - zsi0 ) 
    470          dummy_fac2 = 1._wp / REAL(nlay_i,wp) 
    471  
     468         zfac0 = 1._wp / ( zsi0 - zsi1 ) 
     469         zfac1 = zsi1 / ( zsi1 - zsi0 ) 
    472470         DO jk = 1, nlay_i 
    473471            DO ji = kideb, kiut 
     
    478476               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws  
    479477               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) )  
    480                ! if 2.sm_i GE sss_m then zswibal = 1 
     478               ! if 2.sm_i GE sss_m then rswitch = 1 
    481479               ! this is to force a constant salinity profile in the Baltic Sea 
    482                zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 
     480               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 
    483481               ! 
    484                zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal ) 
     482               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1.0 - rswitch ) 
    485483               ! 
    486                zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2 
     484               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 
    487485               ! weighting the profile 
    488486               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
     
    501499         ! 
    502500         DO jk = 1, nlay_i 
    503             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
     501            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 
    504502            zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  ) 
    505503            DO ji = kideb, kiut 
     
    734732               ! recompute ht_i, ht_s avoiding out of bounds values 
    735733               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 
    736                zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
     734               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 
    737735            ENDIF 
    738736         ENDDO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5076 r5078  
    266266                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 
    267267                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & 
    268                         rswitch / nlay_i 
     268                        rswitch * r1_nlay_i 
    269269                  END DO 
    270270               END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r4990 r5078  
    8282   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
    8383   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
     84#endif 
     85#if defined key_lim3 
     86   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
     87   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
    8488#endif 
    8589   !!---------------------------------------------------------------------- 
     
    166170      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice 
    167171#endif 
    168  
     172#if defined key_lim3 
     173      r1_rhoic = 1._wp / rhoic 
     174      r1_rhosn = 1._wp / rhosn 
     175#endif 
    169176      IF(lwp) THEN 
    170177         WRITE(numout,*) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5070 r5078  
    391391      rdt_ice   = nn_fsbc * rdttra(1)   
    392392      r1_rdtice = 1._wp / rdt_ice  
     393 
     394      ! inverse of nlay_i and nlay_s 
     395      r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 
     396      r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 
    393397      ! 
    394398   END SUBROUTINE ice_run 
Note: See TracChangeset for help on using the changeset viewer.