Changeset 5070


Ignore:
Timestamp:
2015-02-09T14:39:07+01:00 (6 years ago)
Author:
clem
Message:

LIM3 cosmetic changes

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

Legend:

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

    r5067 r5070  
    147147   !! tm_i        |      -      |    Mean sea ice temperature     | K     | 
    148148   !! ot_i        !      -      !    Sea ice areal age content    | day   | 
    149    !! et_i        !      -      !    Total ice enthalpy           | 10^9 J|  
    150    !! et_s        !      -      !    Total snow enthalpy          | 10^9 J|  
     149   !! et_i        !      -      !    Total ice enthalpy           | J/m2  |  
     150   !! et_s        !      -      !    Total snow enthalpy          | J/m2  |  
    151151   !! bv_i        !      -      !    Mean relative brine volume   | ???   |  
    152152   !!===================================================================== 
     
    172172 
    173173   !                                     !!** ice-dynamics namelist (namicedyn) ** 
     174   LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength 
    174175   INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75) 
    175    INTEGER , PUBLIC ::   nn_icestr_bvf    !: use brine volume to diminish ice strength 
    176176   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
     177   INTEGER , PUBLIC ::   nn_ahi0          !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation) 
     178   REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_icestr = 1 
    177179   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
    178180   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength (N/M), Hibler JPO79 
     
    202204   !                                     !!** ice-mechanical redistribution namelist (namiceitdme) 
    203205   REAL(wp), PUBLIC ::   rn_cs            !: fraction of shearing energy contributing to ridging             
    204    REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_ice str = 1 
    205206   REAL(wp), PUBLIC ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging 
    206207   REAL(wp), PUBLIC ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging 
     
    217218 
    218219   !                                     !!** ice-mechanical redistribution namelist (namiceitdme) 
    219    INTEGER , PUBLIC ::   nn_rafting      !: rafting of ice or not                         
     220   LOGICAL , PUBLIC ::   ln_rafting      !: rafting of ice or not                         
    220221   INTEGER , PUBLIC ::   nn_partfun      !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 
    221222 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5067 r5070  
    173173      IF( icount == 0 ) THEN 
    174174 
    175          zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
     175         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     176            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
     177            &                ) *  e12t(:,:) * tmask(:,:,1) ) 
     178 
     179         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     180            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
     181            &                ) *  e12t(:,:) * tmask(:,:,1) ) 
     182 
     183         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     184            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
     185            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
     186 
     187         zvi_b  = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
     188 
    176189         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
    177          zei_b  = glob_sum( SUM( SUM(   e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv +  & 
    178             &               SUM( SUM(   e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
    179          zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    180             &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    181             &             ) * e12t(:,:) * tmask(:,:,1) ) 
    182          zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    183             &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
    184             &                 ) * e12t(:,:) * tmask(:,:,1) ) 
    185          zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    186             &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    187             &                  ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
     190 
     191         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
     192            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
     193                            ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
    188194 
    189195      ELSEIF( icount == 1 ) THEN 
    190196 
    191          zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    192             &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
    193             &                ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 
    194          zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    195             &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    196             &                ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 
    197          zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    198             &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    199             &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 
     197         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     198            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
     199            &              ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 
     200 
     201         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     202            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
     203            &              ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 
     204 
     205         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     206            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
     207            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 
    200208  
    201          zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  & 
     209         zvi  = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  & 
    202210            &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw  
     211 
    203212         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
    204          zei  =   glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv +   & 
    205             &               SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv )  & 
    206             &                         * r1_rdtice - zei_b * r1_rdtice + zft 
    207  
    208          zvmin = glob_min(v_i) 
    209          zamax = glob_max(SUM(a_i,dim=3)) 
    210          zamin = glob_min(a_i) 
     213 
     214         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
     215            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
     216            &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 
     217 
     218         zvmin = glob_min( v_i ) 
     219         zamax = glob_max( SUM( a_i, dim=3 ) ) 
     220         zamin = glob_min( a_i ) 
    211221        
    212222         IF(lwp) THEN 
    213223            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
    214224            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
    215             IF ( ABS( zei    ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
     225            IF ( ABS( zei    ) >  1    ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
    216226            IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
    217227            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

    r5067 r5070  
    180180      DO jj = 1, jpj 
    181181         DO ji = 1, jpi 
    182             IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     182            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    183183               !CALL lim_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    184184               !DO jl = 1, jpl 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r5067 r5070  
    113113 
    114114      ! Heat budget 
    115       zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J] 
    116       zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J] 
     115      zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J] 
     116      zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 
    117117      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    118118      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
     
    137137      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 
    138138      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
    139                           &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 
     139                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 
     140                          &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 
    140141      ! 
    141142      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5067 r5070  
    172172      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    173173         ! 
    174          zcoef = SQRT( 0.5_wp ) / rau0 
     174         zcoef = SQRT( 0.5_wp ) * r1_rau0 
    175175         DO jj = 2, jpjm1 
    176176            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    243243      !!------------------------------------------------------------------- 
    244244      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    245       NAMELIST/namicedyn/ nn_icestr, nn_icestr_bvf, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, nn_nevp, rn_relast, rn_ahi0_ref 
     245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
     246         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 
    246247      INTEGER  ::   ji, jj 
    247248      REAL(wp) ::   za00, zd_max 
     
    261262         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 
    262263         WRITE(numout,*) '~~~~~~~~~~~~' 
    263          WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr =', nn_icestr  
    264          WRITE(numout,*)'    Including brine volume in ice strength comp.     nn_icestr_bvf=', nn_icestr_bvf 
    265          WRITE(numout,*) '   drag coefficient for oceanic stress              rn_cio     = ', rn_cio 
    266          WRITE(numout,*) '   first bulk-rheology parameter                    rn_pstar  = ', rn_pstar 
    267          WRITE(numout,*) '   second bulk-rhelogy parameter                    rn_crhg  = ', rn_crhg 
    268          WRITE(numout,*) '   creep limit                                      rn_creepl = ', rn_creepl 
    269          WRITE(numout,*) '   eccentricity of the elliptical yield curve       rn_ecc    = ', rn_ecc 
    270          WRITE(numout,*) '   number of iterations for subcycling              nn_nevp   = ', nn_nevp 
    271          WRITE(numout,*) '   ratio of elastic timescale over ice time step    rn_relast = ', rn_relast 
    272          WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)       rn_ahi0_ref   = ', rn_ahi0_ref 
     264         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr  
     265         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf 
     266         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging    rn_pe_rdg     = ', rn_pe_rdg  
     267         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio 
     268         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar 
     269         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg 
     270         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl 
     271         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc 
     272         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
     273         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast 
     274         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0 
     275         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref 
    273276      ENDIF 
    274277      ! 
     
    277280      ! 
    278281      !  Diffusion coefficients 
    279       zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    280       IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    281             
    282       za00 = rn_ahi0_ref / ( 1.e05_wp )              ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    283                                               !                    (60° = min latitude for ice cover)   
    284       DO jj = 1, jpj 
    285          DO ji = 1, jpi 
    286             ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
    287             ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     282      SELECT CASE( nn_ahi0 ) 
     283 
     284      CASE( 0 ) 
     285         ahiu(:,:) = rn_ahi0_ref 
     286         ahiv(:,:) = rn_ahi0_ref 
     287 
     288         IF(lwp) WRITE(numout,*) '' 
     289         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
     290 
     291      CASE( 1 )  
     292 
     293         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     294         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
     295          
     296         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
     297                                                        !                    (60° = min latitude for ice cover)   
     298         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
     299 
     300         IF(lwp) WRITE(numout,*) '' 
     301         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
     302         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
     303          
     304      CASE( 2 )  
     305 
     306         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     307         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     308          
     309         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
     310                                                 !                    (60° = min latitude for ice cover)   
     311         DO jj = 1, jpj 
     312            DO ji = 1, jpi 
     313               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
     314               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     315            END DO 
    288316         END DO 
    289       END DO 
    290       ! 
    291       IF(lwp) WRITE(numout,*) '' 
    292       IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
    293       IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
    294   
     317         ! 
     318         IF(lwp) WRITE(numout,*) '' 
     319         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
     320         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
     321          
     322      END SELECT 
     323 
    295324   END SUBROUTINE lim_dyn_init 
    296325 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r5047 r5070  
    7272         DO jj = 2, jpjm1   
    7373            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    74                efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     74               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj) 
    7575            END DO 
    7676         END DO 
     
    9898         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    9999            DO ji = 1 , fs_jpim1   ! vector opt. 
    100                zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    101                zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
     100               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     101               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    102102            END DO 
    103103         END DO 
     
    105105         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    106106            DO ji = fs_2 , fs_jpim1   ! vector opt.  
    107                zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
    108                   &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) ) 
     107               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    109108            END DO 
    110109         END DO 
     
    116115               zrlxint = (   ztab0(ji,jj)    & 
    117116                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   & 
    118                   &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             &  
    119                   &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
     117                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               &  
     118                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
    120119               zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) ) 
    121120            END DO 
     
    139138      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    140139         DO ji = 1 , fs_jpim1   ! vector opt. 
    141             zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
    142             zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
     140            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     141            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
    143142         END DO 
    144143      END DO 
     
    146145      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    147146         DO ji = fs_2 , fs_jpim1   ! vector opt.  
    148             zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
    149                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) ) 
     147            zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 
    150148            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 
    151149         END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5067 r5070  
    157157      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    158158      !-----------------------------------------------------------------------------! 
    159       Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE 
     159      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0                ! proport const for PE 
    160160      ! 
    161161      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
     
    191191            !  (thick, newly ridged ice). 
    192192 
    193             closing_net(ji,jj) = rn_cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     193            closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    194194 
    195195            ! 2.2 divu_adv 
     
    235235               ! Reduce the closing rate if more than 100% of the open water  
    236236               ! would be removed.  Reduce the opening rate proportionately. 
    237                IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
     237               IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN 
    238238                  w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    239                   IF ( w1 .GT. ato_i(ji,jj)) THEN 
     239                  IF ( w1 > ato_i(ji,jj)) THEN 
    240240                     tmpfac = ato_i(ji,jj) / w1 
    241241                     closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
     
    286286         DO jj = 1, jpj 
    287287            DO ji = 1, jpi 
    288                IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN 
     288               IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 
    289289                  closing_net(ji,jj) = 0._wp 
    290290                  opning     (ji,jj) = 0._wp 
     
    357357                  WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
    358358               END DO 
    359             ENDIF  ! asum 
    360  
    361          END DO !ji 
    362       END DO !jj 
     359            ENDIF 
     360         END DO 
     361      END DO 
    363362 
    364363      ! Conservation check 
     
    482481                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )    
    483482!!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                     
    484                   ENDIF            ! aicen > epsi10 
     483                  ENDIF 
    485484                  ! 
    486                END DO ! ji 
    487             END DO !jj 
    488          END DO !jl 
     485               END DO 
     486            END DO 
     487         END DO 
    489488 
    490489         zzc = rn_pe_rdg * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
     
    510509      ! CAN BE REMOVED 
    511510      ! 
    512       IF ( nn_icestr_bvf == 1 ) THEN 
     511      IF( ln_icestr_bvf ) THEN 
    513512 
    514513         DO jj = 1, jpj 
    515514            DO ji = 1, jpi 
    516                IF ( bv_i(ji,jj) .GT. 0.0 ) THEN 
     515               IF ( bv_i(ji,jj) > 0.0 ) THEN 
    517516                  zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 ) 
    518517               ELSE 
     
    520519               ENDIF 
    521520               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 
    522             END DO              ! j 
    523          END DO                 ! i 
     521            END DO 
     522         END DO 
    524523 
    525524      ENDIF 
     
    539538         DO jj = 2, jpj - 1 
    540539            DO ji = 2, jpi - 1 
    541                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 
    542                   ! present 
     540               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 
    543541                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    544542                     &          + strength(ji-1,jj) * tmask(ji-1,jj,1) &   
     
    562560         CALL lbc_lnk( strength, 'T', 1. ) 
    563561 
    564       ENDIF ! ksmooth 
     562      ENDIF 
    565563 
    566564      !-------------------- 
     
    579577         DO jj = 1, jpj - 1 
    580578            DO ji = 1, jpi - 1 
    581                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present 
     579               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN       ! ice is present 
    582580                  numts_rm = 1 ! number of time steps for the running mean 
    583                   IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
    584                   IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
     581                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
     582                  IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
    585583                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 
    586584                  strp2(ji,jj) = strp1(ji,jj) 
     
    661659         DO jj = 1, jpj  
    662660            DO ji = 1, jpi 
    663                IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
    664                ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
     661               IF( a_i(ji,jj,jl) > epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
     662               ELSE                                ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
    665663               ENDIF 
    666664            END DO 
     
    716714      ENDIF ! nn_partfun 
    717715 
    718       IF( nn_rafting == 1 ) THEN      ! Ridging and rafting ice participation functions 
     716      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions 
    719717         ! 
    720718         DO jl = 1, jpl 
    721719            DO jj = 1, jpj  
    722720               DO ji = 1, jpi 
    723                   IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN 
     721                  IF ( athorn(ji,jj,jl) > 0._wp ) THEN 
    724722!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time.... 
    725723                     aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
     
    727725                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp 
    728726                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 
    729                   ENDIF ! athorn 
    730                END DO ! ji 
    731             END DO ! jj 
    732          END DO ! jl 
    733  
    734       ELSE  ! nn_rafting = 0 
     727                  ENDIF 
     728               END DO 
     729            END DO 
     730         END DO 
     731 
     732      ELSE 
    735733         ! 
    736734         DO jl = 1, jpl 
     
    740738      ENDIF 
    741739 
    742       IF ( nn_rafting == 1 ) THEN 
    743  
    744          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 
     740      IF( ln_rafting ) THEN 
     741 
     742         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN 
    745743            DO jl = 1, jpl 
    746744               DO jj = 1, jpj 
    747745                  DO ji = 1, jpi 
    748                      IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN 
     746                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 
    749747                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    750748                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
     
    792790            DO ji = 1, jpi 
    793791 
    794                IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
     792               IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 
    795793                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    796794                  hrmean          = MAX(SQRT(rn_hstar*hi), hi*krdgmin) 
     
    988986         large_afrft = .false. 
    989987 
    990 !CDIR NODEP 
    991988         DO ij = 1, icells 
    992989            ji = indxi(ij) 
     
    11151112         !-------------------------------------------------------------------- 
    11161113         DO jk = 1, nlay_i 
    1117 !CDIR NODEP 
    11181114            DO ij = 1, icells 
    11191115               ji = indxi(ij) 
     
    11421138         IF( con_i ) THEN 
    11431139            DO jk = 1, nlay_i 
    1144 !CDIR NODEP 
    11451140               DO ij = 1, icells 
    11461141                  ji = indxi(ij) 
     
    11521147 
    11531148         IF( large_afrac ) THEN   ! there is a bug 
    1154 !CDIR NODEP 
    11551149            DO ij = 1, icells 
    11561150               ji = indxi(ij) 
     
    11641158         ENDIF 
    11651159         IF( large_afrft ) THEN  ! there is a bug 
    1166 !CDIR NODEP 
    11671160            DO ij = 1, icells 
    11681161               ji = indxi(ij) 
     
    11821175         DO jl2  = 1, jpl  
    11831176            ! over categories to which ridged ice is transferred 
    1184 !CDIR NODEP 
    11851177            DO ij = 1, icells 
    11861178               ji = indxi(ij) 
     
    12151207            ! Transfer ice energy to category jl2 by ridging 
    12161208            DO jk = 1, nlay_i 
    1217 !CDIR NODEP 
    12181209               DO ij = 1, icells 
    12191210                  ji = indxi(ij) 
     
    12271218         DO jl2 = 1, jpl  
    12281219 
    1229 !CDIR NODEP 
    12301220            DO ij = 1, icells 
    12311221               ji = indxi(ij) 
     
    12421232                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
    12431233                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
    1244                ENDIF ! hraft 
     1234               ENDIF 
    12451235               ! 
    1246             END DO ! ij 
     1236            END DO 
    12471237 
    12481238            ! Transfer rafted ice energy to category jl2  
    12491239            DO jk = 1, nlay_i 
    1250 !CDIR NODEP 
    12511240               DO ij = 1, icells 
    12521241                  ji = indxi(ij) 
     
    12561245                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    12571246                  ENDIF 
    1258                END DO           ! ij 
    1259             END DO !jk 
    1260  
    1261          END DO ! jl2 
     1247               END DO 
     1248            END DO 
     1249 
     1250         END DO 
    12621251 
    12631252      END DO ! jl1 (deforming categories) 
     
    13311320      !!------------------------------------------------------------------- 
    13321321      INTEGER :: ios                 ! Local integer output status for namelist read 
    1333       NAMELIST/namiceitdme/ rn_cs, rn_pe_rdg, rn_fsnowrdg, rn_fsnowrft,              &  
    1334         &                   rn_gstar, rn_astar, rn_hstar, nn_rafting, rn_hraft, rn_craft, rn_por_rdg, & 
     1322      NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft,              &  
     1323        &                   rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 
    13351324        &                   nn_partfun 
    13361325      !!------------------------------------------------------------------- 
     
    13491338         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 
    13501339         WRITE(numout,*)' ~~~~~~~~~~~~~~~' 
    1351          WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs              ', rn_cs  
    1352          WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      rn_pe_rdg              ', rn_pe_rdg  
    1353          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg        ', rn_fsnowrdg  
    1354          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft        ', rn_fsnowrft  
    1355          WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar           ', rn_gstar 
    1356          WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar           ', rn_astar 
    1357          WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar           ', rn_hstar 
    1358          WRITE(numout,*)'   Rafting of ice sheets or not                            nn_rafting        ', nn_rafting 
    1359          WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft       ', rn_hraft 
    1360          WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft           ', rn_craft   
    1361          WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg       ', rn_por_rdg 
    1362          WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun     ', nn_partfun 
     1340         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs  
     1341         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg  
     1342         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft  
     1343         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar 
     1344         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar 
     1345         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar 
     1346         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting 
     1347         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft 
     1348         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft   
     1349         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg 
     1350         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun 
    13631351      ENDIF 
    13641352      ! 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5067 r5070  
    160160               zremap_flag(ji,jj) = 0 
    161161            ENDIF 
    162          END DO !ji 
    163       END DO !jj 
     162         END DO 
     163      END DO 
    164164 
    165165      !----------------------------------------------------------------------------------------------- 
     
    201201         END DO 
    202202 
    203       END DO !jl 
     203      END DO 
    204204 
    205205      !----------------------------------------------------------------------------------------------- 
     
    244244      !----------------------------------------------------------------------------------------------- 
    245245      !- 7.1 g(h) for category 1 at start of time step 
    246       CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         & 
    247          &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
     246      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
    248247         &                  hR(:,:,klbnd), zremap_flag ) 
    249248 
     
    254253 
    255254         !ji 
    256          IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 
     255         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
    257256            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    258257            ! ji, a_i > epsi10 
    259             IF (zdh0 .lt. 0.0) THEN !remove area from category 1 
     258            IF( zdh0 < 0.0 ) THEN !remove area from category 1 
    260259               ! ji, a_i > epsi10; zdh0 < 0 
    261                zdh0 = MIN(-zdh0,hi_max(klbnd)) 
     260               zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
    262261 
    263262               !Integrate g(1) from 0 to dh0 to estimate area melted 
    264                zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 
    265                IF (zetamax.gt.0.0) THEN 
     263               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 
     264               IF( zetamax > 0.0 ) THEN 
    266265                  zx1  = zetamax 
    267                   zx2  = 0.5 * zetamax*zetamax  
     266                  zx2  = 0.5 * zetamax * zetamax  
    268267                  zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 
    269268                  ! Constrain new thickness <= ht_i 
    270                   zdamax = a_i(ii,ij,klbnd) * &  
    271                      (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 
     269                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 
    272270                  !ice area lost due to melting of thin ice 
    273                   zda0   = MIN(zda0, zdamax) 
     271                  zda0   = MIN( zda0, zdamax ) 
    274272 
    275273                  ! Remove area, conserving volume 
    276                   ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) &  
    277                      * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
     274                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    278275                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
    279                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 
    280                ENDIF     ! zetamax > 0 
     276                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 
     277               ENDIF 
    281278               ! ji, a_i > epsi10 
    282279 
    283280            ELSE ! if ice accretion 
    284281               ! ji, a_i > epsi10; zdh0 > 0 
    285                zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
     282               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
    286283               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    287284               ! growth in openwater (F0 = f1) 
     
    295292      !- 7.3 g(h) for each thickness category   
    296293      DO jl = klbnd, kubnd 
    297          CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
    298             g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 
     294         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
     295            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
    299296      END DO 
    300297 
     
    316313            ij = nind_j(ji) 
    317314 
    318             IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 
     315            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 
    319316 
    320317               ! left and right integration limits in eta space 
    321                zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 
    322                zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 
     318               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 
     319               zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
    323320               zdonor(ii,ij,jl) = jl 
    324321 
     
    327324               ! left and right integration limits in eta space 
    328325               zvetamin(ji) = 0.0 
    329                zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 
     326               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 
    330327               zdonor(ii,ij,jl) = jl + 1 
    331328 
    332329            ENDIF  ! zhbnew(jl) > hi_max(jl) 
    333330 
    334             zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin 
     331            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 
    335332            zetamin = zvetamin(ji) 
    336333 
    337334            zx1  = zetamax - zetamin 
    338             zwk1 = zetamin*zetamin 
    339             zwk2 = zetamax*zetamax 
    340             zx2  = 0.5 * (zwk2 - zwk1) 
     335            zwk1 = zetamin * zetamin 
     336            zwk2 = zetamax * zetamax 
     337            zx2  = 0.5 * ( zwk2 - zwk1 ) 
    341338            zwk1 = zwk1 * zetamin 
    342339            zwk2 = zwk2 * zetamax 
    343             zx3  = 1.0/3.0 * (zwk2 - zwk1) 
     340            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 ) 
    344341            nd   = zdonor(ii,ij,jl) 
    345342            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
    346343            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    347344 
    348          END DO ! ji 
     345         END DO 
    349346      END DO ! jl klbnd -> kubnd - 1 
    350347 
     
    365362            ht_i(ii,ij,1) = rn_himin 
    366363         ENDIF 
    367       END DO !ji 
     364      END DO 
    368365 
    369366      !!---------------------------------------------------------------------------------------------- 
     
    401398 
    402399 
    403    SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice,   & 
    404       &                        g0, g1, hL, hR, zremap_flag ) 
     400   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
    405401      !!------------------------------------------------------------------ 
    406402      !!                ***  ROUTINE lim_itd_fitline *** 
     
    442438               ! Change hL or hR if hice falls outside central third of range 
    443439 
    444                zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 
    445                zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj)) 
     440               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 
     441               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    446442 
    447443               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) 
     
    454450               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    455451               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 
    456                g0(ji,jj) = zwk1 * ( 2._wp/3._wp - zwk2 ) 
    457                g1(ji,jj) = 2._wp * zdhr * zwk1 * (zwk2 - 0.5) 
     452               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 
     453               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 
    458454               ! 
    459455            ELSE                   ! remap_flag = .false. or a_i < epsi10  
     
    516512 
    517513      DO jl = klbnd, kubnd 
    518          zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 
     514         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 
    519515      END DO 
    520516 
     
    539535            DO ji = 1, jpi 
    540536 
    541                IF (zdonor(ji,jj,jl) .GT. 0) THEN 
     537               IF (zdonor(ji,jj,jl) > 0) THEN 
    542538                  jl1 = zdonor(ji,jj,jl) 
    543539 
    544                   IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 
    545                      IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 
    546                         IF ( ( jl1.EQ.jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )           & 
    547                            .OR.                                      & 
    548                            ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )           &   
    549                            ) THEN                                                              
     540                  IF (zdaice(ji,jj,jl) < 0.0) THEN 
     541                     IF (zdaice(ji,jj,jl) > -epsi10) THEN 
     542                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
     543                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    550544                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
    551545                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
     
    559553                  ENDIF 
    560554 
    561                   IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 
    562                      IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 
    563                         IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) )     & 
    564                            .OR.                                     & 
    565                            ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 
    566                            ) THEN 
     555                  IF (zdvice(ji,jj,jl) < 0.0) THEN 
     556                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 
     557                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
     558                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    567559                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
    568560                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     
    577569 
    578570                  ! If daice is close to aicen, set daice = aicen. 
    579                   IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 
    580                      IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 
     571                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 
     572                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 
    581573                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    582574                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     
    586578                  ENDIF 
    587579 
    588                   IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 
    589                      IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 
     580                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 
     581                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 
    590582                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    591583                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     
    596588 
    597589               ENDIF               ! donor > 0 
    598             END DO                   ! i 
    599          END DO                 ! j 
    600  
    601       END DO !jl 
     590            END DO 
     591         END DO 
     592 
     593      END DO 
    602594 
    603595      !------------------------------------------------------------------------------- 
     
    609601         DO jj = 1, jpj 
    610602            DO ji = 1, jpi 
    611                IF (zdaice(ji,jj,jl) .GT. 0.0 ) THEN ! daice(n) can be < puny 
     603               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 
    612604                  nbrem = nbrem + 1 
    613605                  nind_i(nbrem) = ji 
    614606                  nind_j(nbrem) = jj 
    615                ENDIF ! tmask 
     607               ENDIF 
    616608            END DO 
    617609         END DO 
     
    622614 
    623615            jl1 = zdonor(ii,ij,jl) 
    624             rswitch             = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
    625             zworka(ii,ij)   = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch 
     616            rswitch       = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
     617            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
    626618            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    627             ELSE                    ;   jl2 = jl  
     619            ELSE                  ;   jl2 = jl  
    628620            ENDIF 
    629621 
     
    682674            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
    683675 
    684          END DO                 ! ji 
     676         END DO 
    685677 
    686678         !------------------ 
     
    689681 
    690682         DO jk = 1, nlay_i 
    691 !CDIR NODEP 
    692683            DO ji = 1, nbrem 
    693684               ii = nind_i(ji) 
     
    695686 
    696687               jl1 = zdonor(ii,ij,jl) 
    697                IF (jl1 .EQ. jl) THEN 
     688               IF (jl1 == jl) THEN 
    698689                  jl2 = jl+1 
    699690               ELSE             ! n1 = n+1 
     
    704695               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice 
    705696               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice  
    706             END DO              ! ji 
    707          END DO                 ! jk 
     697            END DO 
     698         END DO 
    708699 
    709700      END DO                   ! boundaries, 1 to ncat-1 
     
    719710                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    720711                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    721                   rswitch         =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
     712                  rswitch         =  1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 
    722713               ELSE 
    723714                  ht_i(ji,jj,jl)  = 0._wp 
    724715                  t_su(ji,jj,jl)  = rtt 
    725716               ENDIF 
    726             END DO                 ! ji 
    727          END DO                 ! jj 
    728       END DO                    ! jl 
     717            END DO 
     718         END DO 
     719      END DO 
    729720      ! 
    730721      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 
     
    836827                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 
    837828               ENDIF 
    838             END DO                 ! ji 
    839          END DO                 ! jj 
     829            END DO 
     830         END DO 
    840831         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
    841832 
     
    870861                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
    871862               ENDIF 
    872             END DO                 ! ji 
    873          END DO                 ! jj 
     863            END DO 
     864         END DO 
    874865 
    875866         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     
    890881!                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    891882!               ENDIF 
    892 !            END DO                 ! ji 
    893 !         END DO                 ! jj 
     883!            END DO 
     884!         END DO 
    894885         ! clem-change end 
    895886 
    896       END DO                    ! jl 
     887      END DO 
    897888 
    898889      !------------------------------------------------------------------------------ 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5067 r5070  
    291291            ! include it later 
    292292 
    293             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    294             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     293            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     294            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    295295 
    296296            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
     
    357357                  &            ) * r1_e12t(ji,jj) 
    358358 
    359                zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    360                   &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     359               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     360                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    361361                  &         ) * r1_e12t(ji,jj) 
    362362 
    363363               ! 
    364                zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    365                   &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     364               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     365                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    366366                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   & 
    367367                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
     
    391391 
    392392               !- Calculate Delta on corners 
    393                zddc  = (  ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    394                   &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     393               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     394                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    395395                  &    ) * r1_e12f(ji,jj) 
    396396 
    397                zdtc  = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    398                   &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     397               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     398                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    399399                  &    ) * r1_e12f(ji,jj) 
    400400 
     
    420420               !- contribution of zs1, zs2 and zs12 to zf1 
    421421               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  & 
    422                   &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj)          & 
    423                   &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj)  & 
     422                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          & 
     423                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  & 
    424424                  &                ) * r1_e12u(ji,jj) 
    425425               ! contribution of zs1, zs2 and zs12 to zf2 
    426426               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
    427                   &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj)          & 
    428                   &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj)  & 
     427                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          & 
     428                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  & 
    429429                  &               )  * r1_e12v(ji,jj) 
    430430            END DO 
     
    609609                  &            ) * r1_e12t(ji,jj) 
    610610 
    611                zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
    612                   &          -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
     611               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
     612                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
    613613                  &         ) * r1_e12t(ji,jj) 
    614614               ! 
    615615               ! SB modif because ocean has no slip boundary condition  
    616                zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    617                   &          +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     616               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     617                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    618618                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     & 
    619619                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
     
    677677            DO jj = k_j1+1, k_jpj-1 
    678678               DO ji = 2, jpim1 
    679                   IF (zpresh(ji,jj) .GT. 1.0) THEN 
     679                  IF (zpresh(ji,jj) > 1.0) THEN 
    680680                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    681681                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5067 r5070  
    631631            rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 
    632632            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 
    633             ! adjust thickness 
     633            ! retrieve total concentration 
    634634            at_i_1d(ji) = a_i_1d(ji) 
    635635         END IF 
     
    651651      !!------------------------------------------------------------------- 
    652652      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    653       NAMELIST/namicethd/ rn_hnewice, nn_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,    & 
    654          &                rn_himin, parsub, rn_betas,                          &  
    655          &                rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
     653      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
     654         &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
    656655         &                nn_monocat 
    657656      !!------------------------------------------------------------------- 
     
    682681         WRITE(numout,*) 
    683682         WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation ' 
    684          WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice      = ', rn_hnewice 
    685          WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       nn_frazil     = ', nn_frazil 
    686          WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb     = ', rn_maxfrazb 
    687          WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb       = ', rn_vfrazb 
    688          WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb       = ', rn_Cfrazb 
    689          WRITE(numout,*)'      minimum ice thickness                                   rn_himin       = ', rn_himin  
     683         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice 
     684         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil 
     685         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb 
     686         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb    = ', rn_vfrazb 
     687         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb 
     688         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin  
    690689         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice ' 
    691690         WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub   
    692          WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas        = ', rn_betas 
    693          WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i      = ', rn_kappa_i 
     691         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas 
     692         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
    694693         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nn_conv_dif  = ', nn_conv_dif 
    695694         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif 
    696          WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon  = ', nn_ice_thcon 
     695         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon 
    697696         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i 
    698697         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5067 r5070  
    276276 
    277277            ! updates available heat + thickness 
    278             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
     278            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    279279            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    280280 
     
    573573            
    574574            ENDIF 
    575          END DO ! ji 
    576       END DO ! jk 
     575         END DO 
     576      END DO 
    577577 
    578578      !------------------------------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5067 r5070  
    517517 
    518518         DO ji = kideb , kiut 
    519             IF ( ht_s_1d(ji).gt.0.0 ) THEN 
     519            IF ( ht_s_1d(ji) > 0.0 ) THEN 
    520520               ! 
    521521               !------------------------------------------------------------------------------| 
     
    541541               ENDIF 
    542542 
    543                IF ( t_su_1d(ji) .LT. rtt ) THEN 
     543               IF ( t_su_1d(ji) < rtt ) THEN 
    544544 
    545545                  !------------------------------------------------------------------------------| 
     
    587587               !------------------------------------------------------------------------------| 
    588588               ! 
    589                IF (t_su_1d(ji) .LT. rtt) THEN 
     589               IF (t_su_1d(ji) < rtt) THEN 
    590590                  ! 
    591591                  !------------------------------------------------------------------------------| 
     
    705705         DO ji = kideb , kiut 
    706706            ! snow temperatures       
    707             IF (ht_s_1d(ji).GT.0._wp) & 
     707            IF (ht_s_1d(ji) > 0._wp) & 
    708708               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    709709               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5067 r5070  
    156156      hicol(:,:) = rn_hnewice 
    157157 
    158       IF( nn_frazil == 1 ) THEN 
     158      IF( ln_frazil ) THEN 
    159159 
    160160         !-------------------- 
     
    221221                  iterate_frazil = .true. 
    222222 
    223                   DO WHILE ( iter .LT. 100 .AND. iterate_frazil )  
     223                  DO WHILE ( iter < 100 .AND. iterate_frazil )  
    224224                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 
    225225                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 
     
    321321            zh_newice(ji) = rn_hnewice 
    322322         END DO 
    323          IF( nn_frazil == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
     323         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
    324324 
    325325         !---------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r5067 r5070  
    150150         WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity ' 
    151151         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    152          WRITE(numout,*) ' switch for salinity nn_icesal        : ', nn_icesal 
    153          WRITE(numout,*) ' bulk salinity value if nn_icesal = 1 : ', rn_icesal 
    154          WRITE(numout,*) ' restoring salinity for GD          : ', rn_sal_gd 
    155          WRITE(numout,*) ' restoring time for GD              : ', rn_time_gd 
    156          WRITE(numout,*) ' restoring salinity for flushing    : ', rn_sal_fl 
    157          WRITE(numout,*) ' restoring time for flushing        : ', rn_time_fl 
    158          WRITE(numout,*) ' Maximum tolerated ice salinity     : ', rn_simax 
    159          WRITE(numout,*) ' Minimum tolerated ice salinity     : ', rn_simin 
     152         WRITE(numout,*) '   switch for salinity nn_icesal        = ', nn_icesal 
     153         WRITE(numout,*) '   bulk salinity value if nn_icesal = 1 = ', rn_icesal 
     154         WRITE(numout,*) '   restoring salinity for GD            = ', rn_sal_gd 
     155         WRITE(numout,*) '   restoring time for GD                = ', rn_time_gd 
     156         WRITE(numout,*) '   restoring salinity for flushing      = ', rn_sal_fl 
     157         WRITE(numout,*) '   restoring time for flushing          = ', rn_time_fl 
     158         WRITE(numout,*) '   Maximum tolerated ice salinity       = ', rn_simax 
     159         WRITE(numout,*) '   Minimum tolerated ice salinity       = ', rn_simin 
    160160      ENDIF 
    161161      ! 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5067 r5070  
    106106         ! mass and salt flux init 
    107107         zviold(:,:,:) = v_i(:,:,:) 
     108         zvsold(:,:,:) = v_s(:,:,:) 
    108109         zeiold(:,:)   = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
    109110         zesold(:,:)   = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
     
    131132 
    132133         ! If ice drift field is too fast, use an appropriate time step for advection.          
    133          zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )         ! CFL test for stability 
    134          zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 
     134         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability 
     135         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    135136         IF(lk_mpp )   CALL mpp_max( zcfl ) 
    136137 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5067 r5070  
    2626   REAL(wp), PUBLIC ::   rn_hnewice  !: thickness for new ice formation (m) 
    2727 
    28    INTEGER , PUBLIC ::   nn_frazil   !: use of frazil ice collection as function of wind (1) or not (0) 
     28   LOGICAL , PUBLIC ::   ln_frazil   !: use of frazil ice collection as function of wind (T) or not (F) 
    2929 
    3030   !!----------------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5067 r5070  
    369369         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
    370370         WRITE(numout,*) ' ~~~~~~' 
    371          WRITE(numout,*) '   number of ice  categories                                               = ', jpl 
    372          WRITE(numout,*) '   number of ice  layers                                                   = ', nlay_i 
    373          WRITE(numout,*) '   number of snow layers                                                   = ', nlay_s 
     371         WRITE(numout,*) '   number of ice  categories                               = ', jpl 
     372         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i 
     373         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
    374374         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    375375         WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax  
Note: See TracChangeset for help on using the changeset viewer.