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

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

Ignore:
Timestamp:
2015-01-20T15:26:13+01:00 (9 years ago)
Author:
jamesharle
Message:

Merging branch with HEAD of the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4346 r5038  
    55   !!====================================================================== 
    66   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code  
    7    !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 
     7   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn 
    88   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
     
    2222   USE limthd_lac       ! LIM 
    2323   USE limvar           ! LIM 
    24    USE limcons          ! LIM 
    2524   USE in_out_manager   ! I/O manager 
    2625   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     
    3029  ! Check budget (Rousset) 
    3130   USE iom              ! I/O manager 
    32    USE lib_fortran     ! glob_sum 
     31   USE lib_fortran      ! glob_sum 
    3332   USE limdiahsb 
    34    USE timing          ! Timing 
     33   USE timing           ! Timing 
     34   USE limcons          ! conservation tests 
    3535 
    3636   IMPLICIT NONE 
     
    4242   PUBLIC   lim_itd_me_zapsmall 
    4343   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
    44  
    45    REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    46    REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values 
    47    REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values 
    4844 
    4945   !----------------------------------------------------------------------- 
     
    143139      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    144140      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    145       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
    146       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    147       ! mass and salt flux (clem) 
    148       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
     141      ! 
     142      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    149143      !!----------------------------------------------------------------------------- 
    150144      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    151145 
    152146      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    153  
    154       CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    155  
    156       IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only) 
    157147 
    158148      IF(ln_ctl) THEN 
     
    162152 
    163153      IF( ln_limdyn ) THEN          !   Start ridging and rafting   ! 
    164       ! ------------------------------- 
    165       !- check conservation (C Rousset) 
    166       IF (ln_limdiahsb) THEN 
    167          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    168          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    169          zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
    170          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
    171       ENDIF 
    172       !- check conservation (C Rousset) 
    173       ! ------------------------------- 
    174  
    175       ! mass and salt flux init (clem) 
    176       zviold(:,:,:) = v_i(:,:,:) 
    177       zvsold(:,:,:) = v_s(:,:,:) 
    178       zsmvold(:,:,:) = smv_i(:,:,:) 
     154 
     155      ! conservation test 
     156      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    179157 
    180158      !-----------------------------------------------------------------------------! 
     
    362340            ! 5) Heat, salt and freshwater fluxes 
    363341            !-----------------------------------------------------------------------------! 
    364             fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    365             fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean 
     342            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     343            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2) 
    366344 
    367345         END DO 
     
    399377      CALL lim_itd_me_zapsmall 
    400378 
    401       !-------------------------------- 
    402       ! Update mass/salt fluxes (clem) 
    403       !-------------------------------- 
    404       DO jl = 1, jpl 
    405          DO jj = 1, jpj  
    406             DO ji = 1, jpi 
    407                diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 
    408                rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic  
    409                rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn  
    410                sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice  
    411             END DO 
    412          END DO 
    413       END DO 
    414379 
    415380      IF(ln_ctl) THEN     ! Control print 
     
    445410      ENDIF 
    446411 
    447       ! ------------------------------- 
    448       !- check conservation (C Rousset) 
    449       IF (ln_limdiahsb) THEN 
    450          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    451          zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    452   
    453          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
    454          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    455  
    456          zchk_vmin = glob_min(v_i) 
    457          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    458          zchk_amin = glob_min(a_i) 
    459         
    460          IF(lwp) THEN 
    461             IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday) 
    462             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 
    463             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
    464             IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
    465             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin 
    466          ENDIF 
    467       ENDIF 
    468       !- check conservation (C Rousset) 
    469       ! ------------------------------- 
     412      ! conservation test 
     413      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    470414 
    471415      ENDIF  ! ln_limdyn=.true. 
    472416      ! 
    473417      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    474       ! 
    475       CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    476418      ! 
    477419      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
     
    670612      !!---------------------------------------------------------------------! 
    671613      INTEGER ::   ji,jj, jl    ! dummy loop indices 
    672       INTEGER ::   krdg_index   !  
    673614      REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar 
    674615      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here 
     
    746687      !----------------------------------------------------------------- 
    747688 
    748       krdg_index = 1 
    749  
    750       IF( krdg_index == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
    751          DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates 
     689      IF( partfun_swi == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
     690         DO jl = 0, jpl     
    752691            DO jj = 1, jpj  
    753692               DO ji = 1, jpi 
     
    772711            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
    773712         END DO !jl 
    774          DO jl = 0, ice_cat_bounds(1,2) 
     713         DO jl = 0, jpl 
    775714             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
    776715         END DO 
    777716         ! 
    778       ENDIF ! krdg_index 
    779  
    780       IF( raftswi == 1 ) THEN      ! Ridging and rafting ice participation functions 
     717      ENDIF ! partfun_swi 
     718 
     719      IF( raft_swi == 1 ) THEN      ! Ridging and rafting ice participation functions 
    781720         ! 
    782721         DO jl = 1, jpl 
     
    794733         END DO ! jl 
    795734 
    796       ELSE  ! raftswi = 0 
     735      ELSE  ! raft_swi = 0 
    797736         ! 
    798737         DO jl = 1, jpl 
     
    802741      ENDIF 
    803742 
    804       IF ( raftswi == 1 ) THEN 
     743      IF ( raft_swi == 1 ) THEN 
    805744 
    806745         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 
     
    908847      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    909848      INTEGER ::   icells            ! number of cells with aicen > puny 
    910       REAL(wp) ::   zindb, zsrdg2   ! local scalar 
    911849      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
     850      REAL(wp) ::   zsstK            ! SST in Kelvin 
    912851 
    913852      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
     
    917856 
    918857      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging 
    919       REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging 
     858      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging 
    920859      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging 
    921860 
     
    952891      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    953892      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    954       CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init ) 
    955       CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw ) 
    956       CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init ) 
     893      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
     894      CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
     895      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
    957896 
    958897      ! Conservation check 
     
    1008947         aicen_init(:,:,jl) = a_i(:,:,jl) 
    1009948         vicen_init(:,:,jl) = v_i(:,:,jl) 
    1010          vsnon_init(:,:,jl) = v_s(:,:,jl) 
     949         vsnwn_init(:,:,jl) = v_s(:,:,jl) 
    1011950         ! 
    1012951         smv_i_init(:,:,jl) = smv_i(:,:,jl) 
     
    1014953      END DO !jl 
    1015954 
    1016       esnon_init(:,:,:) = e_s(:,:,1,:) 
     955      esnwn_init(:,:,:) = e_s(:,:,1,:) 
    1017956 
    1018957      DO jl = 1, jpl   
     
    10911030            !     / rafting category n1. 
    10921031            !-------------------------------------------------------------------------- 
    1093             vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
     1032            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 
    10941033            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 
    10951034            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por 
    10961035 
    1097             vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj) 
    1098             esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 
    1099             srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
    1100             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     1036            vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     1037            esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     1038            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     1039            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 
    11011040 
    11021041            ! rafting volumes, heat contents ... 
    11031042            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    1104             vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj) 
    1105             esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) 
     1043            vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     1044            esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    11061045            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
    11071046 
     
    11201059            ! Salinity 
    11211060            !------------- 
    1122             smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids 
    1123  
    1124             zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    1125  
    1126             srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
     1061            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
     1062            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
     1063 
     1064            !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
    11271065             
    1128             !                                                             ! excess of salt is flushed into the ocean 
    1129             !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
    1130  
    1131             !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids              
     1066            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
     1067            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
    11321068 
    11331069            !------------------------------------             
     
    11581094               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    11591095 
    1160             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
    1161                &                                + esrft(ji,jj)*(1.0-fsnowrft)           
     1096            ! in 1e-9 Joules (same as e_s) 
     1097            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
     1098               &                                - esrft(ji,jj)*(1.0-fsnowrft)           
    11621099 
    11631100            !----------------------------------------------------------------- 
     
    11841121               jj = indxj(ij) 
    11851122               ! heat content of ridged ice 
    1186                erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )  
     1123               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)  
    11871124               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
    11881125               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 
    1189                ! sea water heat content 
    1190                ztmelts          = - tmut * sss_m(ji,jj) + rtt 
    1191                ! heat content per unit volume 
    1192                zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 
    1193  
    1194                ! corrected sea water salinity 
    1195                zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 
    1196                zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 
    1197  
    1198                ztmelts          = - tmut * zdummy + rtt 
    1199                ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 
    1200  
    1201                ! heat flux 
    1202                fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 
     1126                
     1127                
     1128               ! enthalpy of the trapped seawater (J/m2, >0) 
     1129               ! clem: if sst>0, then ersw <0 (is that possible?) 
     1130               zsstK  = sst_m(ji,jj) + rt0 
     1131               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 
     1132 
     1133               ! heat flux to the ocean 
     1134               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    12031135 
    12041136               ! Correct dimensions to avoid big values 
    1205                ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09 
    1206  
    1207                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    1208                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 
     1137               ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac 
     1138 
     1139               ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 
     1140               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean  
     1141               !! MV HC 2014 
     1142               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) 
    12091143 
    12101144               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     1145 
    12111146            END DO ! ij 
    12121147         END DO !jk 
     
    12531188         !------------------------------------------------------------------------------- 
    12541189         !        jl1 looping 1-jpl 
    1255          DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2)  
     1190         DO jl2  = 1, jpl  
    12561191            ! over categories to which ridged ice is transferred 
    12571192!CDIR NODEP 
     
    12981233         END DO                 ! jl2 (new ridges)             
    12991234 
    1300          DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)  
     1235         DO jl2 = 1, jpl  
    13011236 
    13021237!CDIR NODEP 
     
    13611296      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    13621297      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    1363       CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init ) 
    1364       CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw ) 
    1365       CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init ) 
     1298      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
     1299      CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
     1300      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
    13661301      ! 
    13671302   END SUBROUTINE lim_itd_me_ridgeshift 
     
    14041339      !!------------------------------------------------------------------- 
    14051340      INTEGER :: ios                 ! Local integer output status for namelist read 
    1406       NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,&  
    1407          Gstar, astar,                                & 
    1408          Hstar, raftswi, hparmeter, Craft, ridge_por, & 
    1409          sal_max_ridge,  partfun_swi, transfun_swi,   & 
    1410          brinstren_swi 
     1341      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,              &  
     1342        &                   Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, & 
     1343        &                   partfun_swi, brinstren_swi 
    14111344      !!------------------------------------------------------------------- 
    14121345      ! 
     
    14181351      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 ) 
    14191352902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp ) 
    1420       WRITE ( numoni, namiceitdme ) 
     1353      IF(lwm) WRITE ( numoni, namiceitdme ) 
    14211354      ! 
    14221355      IF (lwp) THEN                          ! control print 
     
    14321365         WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar 
    14331366         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar 
    1434          WRITE(numout,*)'   Rafting of ice sheets or not                            raftswi         ', raftswi 
     1367         WRITE(numout,*)'   Rafting of ice sheets or not                            raft_swi        ', raft_swi 
    14351368         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter 
    14361369         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft   
    14371370         WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por 
    1438          WRITE(numout,*)'   Maximum salinity of ridging ice                         sal_max_ridge   ', sal_max_ridge 
    14391371         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi 
    1440          WRITE(numout,*)'   Switch for tran. function (0) linear (1) exponential    transfun_swi    ', transfun_swi 
    14411372         WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi 
    14421373      ENDIF 
     
    14621393 
    14631394      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1464       REAL(wp)                          ::   zmask_glo 
     1395      REAL(wp)                          ::   zmask_glo, zsal, zvi, zvs, zei, zes 
    14651396!!gm      REAL(wp) ::   xtmp      ! temporary variable 
    14661397      !!------------------------------------------------------------------- 
     
    14681399      CALL wrk_alloc( jpi, jpj, zmask ) 
    14691400 
     1401      ! to be sure that at_i is the sum of a_i(jl) 
     1402      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     1403 
    14701404      DO jl = 1, jpl 
    1471  
    14721405         !----------------------------------------------------------------- 
    14731406         ! Count categories to be zapped. 
    1474          ! Abort model in case of negative area. 
    14751407         !----------------------------------------------------------------- 
    14761408         icells = 0 
     
    14781410         DO jj = 1, jpj 
    14791411            DO ji = 1, jpi 
    1480                IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
    1481                   & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
    1482                   & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
    1483                   & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1412               IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 
     1413                  zmask(ji,jj) = 1._wp 
     1414               ENDIF 
    14841415            END DO 
    14851416         END DO 
     
    14941425            DO jj = 1 , jpj 
    14951426               DO ji = 1 , jpi 
    1496 !!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice 
    1497 !!gm                  xtmp = xtmp * unit_fac 
    1498                   ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
     1427                  zei  = e_i(ji,jj,jk,jl) 
    14991428                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
     1429                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 
     1430                  ! update exchanges with ocean 
     1431                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    15001432               END DO 
    15011433            END DO 
     
    15041436         DO jj = 1 , jpj 
    15051437            DO ji = 1 , jpi 
    1506  
     1438                
     1439               zsal = smv_i(ji,jj,jl) 
     1440               zvi  = v_i(ji,jj,jl) 
     1441               zvs  = v_s(ji,jj,jl) 
     1442               zes  = e_s(ji,jj,1,jl) 
    15071443               !----------------------------------------------------------------- 
    15081444               ! Zap snow energy and use ocean heat to melt snow 
     
    15141450               ! fluxes are positive to the ocean 
    15151451               ! here the flux has to be negative for the ocean 
    1516 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice 
    1517                !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1518  
    1519 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    1520  
    15211452               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    15221453 
     
    15241455               ! zap ice and snow volume, add water and salt to ocean 
    15251456               !----------------------------------------------------------------- 
    1526  
    1527                !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 
    1528                !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   & 
    1529                !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice 
    1530                !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &  
    1531                !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice 
    1532                !           sfx (i,j)      = sfx (i,j)      + xtmp 
    1533  
    1534                ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj) 
     1457               ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj) 
    15351458               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    15361459               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     
    15391462               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    15401463               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1541                ! 
     1464               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
     1465               ! additional condition 
     1466               IF( v_s(ji,jj,jl) <= epsi10 ) THEN 
     1467                  v_s(ji,jj,jl)   = 0._wp 
     1468                  e_s(ji,jj,1,jl) = 0._wp 
     1469               ENDIF 
     1470               ! update exchanges with ocean 
     1471               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     1472               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
     1473               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
     1474               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    15421475            END DO 
    15431476         END DO 
    1544          ! 
    1545       END DO                 ! jl  
     1477      END DO ! jl  
     1478 
     1479      ! to be sure that at_i is the sum of a_i(jl) 
     1480      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    15461481      ! 
    15471482      CALL wrk_dealloc( jpi, jpj, zmask ) 
Note: See TracChangeset for help on using the changeset viewer.