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 3517 for branches/2012/dev_r3385_NOCS04_HAMF – NEMO

Ignore:
Timestamp:
2012-10-26T12:13:21+02:00 (12 years ago)
Author:
gm
Message:

gm: Branch: dev_r3385_NOCS04_HAMF; #665. update sbccpl ; change LIM3 from equivalent salt flux to salt flux and mass flux

Location:
branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90

    r3489 r3517  
    199199            ! 
    200200            !                          !------------------------------------------! 
    201             !                          !      mass flux at the ocean surface      ! 
     201            !                          !  mass and salt flux at the ocean surface ! 
    202202            !                          !------------------------------------------! 
    203203            ! 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r3396 r3517  
    158158   !! * Share Module variables 
    159159   !!-------------------------------------------------------------------------- 
    160    INTEGER , PUBLIC ::   nstart    !: iteration number of the begining of the run  
    161    INTEGER , PUBLIC ::   nlast     !: iteration number of the end of the run  
    162    INTEGER , PUBLIC ::   nitrun    !: number of iteration 
    163    INTEGER , PUBLIC ::   numit     !: iteration number 
    164    REAL(wp), PUBLIC ::   rdt_ice   !: ice time step 
     160   INTEGER , PUBLIC ::   nstart      !: iteration number of the begining of the run  
     161   INTEGER , PUBLIC ::   nlast       !: iteration number of the end of the run  
     162   INTEGER , PUBLIC ::   nitrun      !: number of iteration 
     163   INTEGER , PUBLIC ::   numit       !: iteration number 
     164   REAL(wp), PUBLIC ::   rdt_ice     !: ice time step 
     165   REAL(wp), PUBLIC ::   r1_rdtice   !: = 1. / rdt_ice 
    165166 
    166167   !                                          !!** ice-dynamic namelist (namicedyn) ** 
    … …  
    201202   !                                              !!** ice-salinity namelist (namicesal) ** 
    202203   INTEGER , PUBLIC ::   num_sal     = 1           !: salinity configuration used in the model 
    203    !                                               ! 1 - s constant in space and time 
     204   !                                               ! 1 - constant salinity in both space and time 
    204205   !                                               ! 2 - prognostic salinity (s(z,t)) 
    205206   !                                               ! 3 - salinity profile, constant in time 
    206    !                                               ! 4 - salinity variations affect only ice thermodynamics 
    207207   INTEGER , PUBLIC ::   sal_prof    = 1           !: salinity profile or not  
    208208   INTEGER , PUBLIC ::   thcon_i_swi = 1           !: thermal conductivity: =1 Untersteiner (1964) ; =2 Pringle et al (2007) 
    … …  
    278278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qfvbq       !: store energy in case of total lateral ablation (?) 
    279279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dmgwi       !: Variation of the mass of snow ice 
    280    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness 
    281    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection 
    282    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation 
    283    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation 
     280   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness [???] 
     281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection                      [???] 
     282   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation [???] 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation [???] 
    284284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection 
    285285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: Mass flux due to snow loss during compression 
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: Equivalent salt flux due to ice growth/melt 
     286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: salt flux due to ice growth/melt   [units ???] 
    287287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhmec       !: Heat flux due to snow loss during compression 
    288288   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: Residual heat flux due to correction of ice thickness 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r3294 r3517  
    7979      CALL lim_thd_sal_init            ! set ice salinity parameters 
    8080      ! 
    81       rdt_ice = nn_fsbc * rdttra(1)    ! sea-ice timestep 
     81      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep 
     82      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse 
    8283      ! 
    8384      CALL lim_msh                     ! ice mesh initialization 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r3294 r3517  
    8888            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    8989               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
    90             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     90            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    9191 
    9292            ps0 (ji,jj) = zslpmax   
    … …  
    273273            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    274274               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    275             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     275            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    276276            ! 
    277277            ps0 (ji,jj) = zslpmax   

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90

    r3402 r3517  
    180180               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp  
    181181               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp 
    182                vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW 
     182               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice  ! volume acc in OW 
    183183            ENDIF 
    184184         END DO 
    … …  
    240240      vinfor(84) = 0.0 
    241241      DO ji = 134, 138 
    242          vinfor(83) = vinfor(83) - v_ice(ji,jj) * &  
    243             e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 
    244          vinfor(84) = vinfor(84) - v_ice(ji,jj) * &  
    245             e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 
     242         vinfor(83) = vinfor(83) - v_ice(ji,jj) * e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 
     243         vinfor(84) = vinfor(84) - v_ice(ji,jj) * e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 
    246244      END DO 
    247245 
    … …  
    331329               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp  
    332330               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp 
    333                vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW 
     331               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice  ! volume acc in OW 
    334332            ENDIF 
    335333         END DO 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r3402 r3517  
    3838   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
    3939 
    40    REAL(wp)  ::   epsi11 = 1.e-11_wp   ! constant values 
    41    REAL(wp)  ::   epsi10 = 1.e-10_wp   ! constant values 
    42    REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
     40   REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values 
     41   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values 
     42   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values 
    4343 
    4444   !----------------------------------------------------------------------- 
    … …  
    4747   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    4848   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    49  
     49   ! 
    5050   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
    5151   !                                                           !  closing associated w/ category n 
    52  
     52   ! 
    5353   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    5454   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    … …  
    7070   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s) 
    7171   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s) 
     72   ! 
    7273   !!---------------------------------------------------------------------- 
    7374   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    … …  
    126127      INTEGER ::   ji, jj, jk, jl   ! dummy loop index 
    127128      INTEGER ::   niter, nitermax = 20   ! local integer  
    128       LOGICAL  ::   asum_error              ! flag for asum .ne. 1 
    129       INTEGER  ::   iterate_ridging         ! if true, repeat the ridging 
    130       REAL(wp) ::   w1, tmpfac, dti         ! local scalar 
     129      LOGICAL  ::   asum_error            ! flag for asum .ne. 1 
     130      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging 
     131      REAL(wp) ::   w1, tmpfac            ! local scalar 
    131132      CHARACTER (len = 15) ::   fieldid 
    132133      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    … …  
    152153      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    153154      !-----------------------------------------------------------------------------! 
    154       ! Set hi_max(ncat) to a big value to ensure that all ridged ice  
    155       ! is thinner than hi_max(ncat). 
     155      ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat). 
    156156 
    157157      hi_max(jpl) = 999.99 
    158158 
    159       Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0      ! proport const for PE 
    160       CALL lim_itd_me_ridgeprep ! prepare ridging 
    161  
     159      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE 
     160      ! 
     161      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
     162      ! 
    162163      IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check 
    163164 
    … …  
    166167            msnow_mlt(ji,jj) = 0._wp 
    167168            esnow_mlt(ji,jj) = 0._wp 
    168             dardg1dt (ji,jj)  = 0._wp 
    169             dardg2dt (ji,jj)  = 0._wp 
    170             dvirdgdt (ji,jj)  = 0._wp 
    171             opening  (ji,jj)  = 0._wp 
     169            dardg1dt (ji,jj) = 0._wp 
     170            dardg2dt (ji,jj) = 0._wp 
     171            dvirdgdt (ji,jj) = 0._wp 
     172            opening  (ji,jj) = 0._wp 
    172173 
    173174            !-----------------------------------------------------------------------------! 
    … …  
    201202            ! to give asum = 1.0 after ridging. 
    202203 
    203             divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) / rdt_ice  ! asum found in ridgeprep 
     204            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    204205 
    205206            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
    … …  
    207208            ! 2.3 opning 
    208209            !------------ 
    209             ! Compute the (non-negative) opening rate that will give  
    210             ! asum = 1.0 after ridging. 
     210            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 
    211211            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 
    212212         END DO 
    … …  
    257257                  IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    258258                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    259                      IF ( w1 > a_i(ji,jj,jl) ) THEN 
     259                     IF ( w1  > a_i(ji,jj,jl) ) THEN 
    260260                        tmpfac = a_i(ji,jj,jl) / w1 
    261261                        closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    … …  
    291291               ELSE 
    292292                  iterate_ridging    = 1 
    293                   divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) / rdt_ice 
     293                  divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) * r1_rdtice 
    294294                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    295295                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
    … …  
    308308 
    309309         IF( iterate_ridging == 1 ) THEN 
    310             IF( niter .GT. nitermax ) THEN 
     310            IF( niter  > nitermax ) THEN 
    311311               WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
    312312               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
    … …  
    323323      ! Update fresh water and heat fluxes due to snow melt. 
    324324 
    325       dti = 1._wp / rdt_ice 
    326  
    327325      asum_error = .false.  
    328326 
    … …  
    330328         DO ji = 1, jpi 
    331329 
    332             IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 
    333  
    334             dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 
    335             dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 
    336             dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 
    337             opening (ji,jj) = opening (ji,jj) * dti 
     330            IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )  asum_error = .true. 
     331 
     332            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     333            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice 
     334            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice 
     335            opening (ji,jj) = opening (ji,jj) * r1_rdtice 
    338336 
    339337            !-----------------------------------------------------------------------------! 
    340338            ! 5) Heat, salt and freshwater fluxes 
    341339            !-----------------------------------------------------------------------------! 
    342             fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * dti     ! fresh water source for ocean 
    343             fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * dti     ! heat sink for ocean 
     340            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     341            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean 
    344342 
    345343         END DO 
    … …  
    349347      DO jj = 1, jpj 
    350348         DO ji = 1, jpi 
    351             IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug 
     349            IF( ABS( asum(ji,jj) - 1._wp )  >  epsi11 ) THEN  ! there is a bug 
    352350               WRITE(numout,*) ' ' 
    353351               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    … …  
    391389      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    392390      d_smv_i_trp(:,:,:)   = 0._wp 
    393       IF(  num_sal == 2  .OR.  num_sal == 4  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
     391      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    394392 
    395393      IF(ln_ctl) THEN     ! Control print 
    … …  
    430428 
    431429      ! update of fields will be made later in lim update 
    432       u_ice(:,:)    = old_u_ice(:,:) 
    433       v_ice(:,:)    = old_v_ice(:,:) 
    434       a_i(:,:,:)    = old_a_i(:,:,:) 
    435       v_s(:,:,:)    = old_v_s(:,:,:) 
    436       v_i(:,:,:)    = old_v_i(:,:,:) 
    437       e_s(:,:,:,:)  = old_e_s(:,:,:,:) 
    438       e_i(:,:,:,:)  = old_e_i(:,:,:,:) 
    439       oa_i(:,:,:)   = old_oa_i(:,:,:) 
    440       IF(  num_sal == 2  .OR.  num_sal == 4  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
     430      u_ice(:,:)   = old_u_ice(:,:) 
     431      v_ice(:,:)   = old_v_ice(:,:) 
     432      a_i(:,:,:)   = old_a_i(:,:,:) 
     433      v_s(:,:,:)   = old_v_s(:,:,:) 
     434      v_i(:,:,:)   = old_v_i(:,:,:) 
     435      e_s(:,:,:,:) = old_e_s(:,:,:,:) 
     436      e_i(:,:,:,:) = old_e_i(:,:,:,:) 
     437      oa_i(:,:,:)  = old_oa_i(:,:,:) 
     438      IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    441439 
    442440      !----------------------------------------------------! 
    … …  
    465463         DO jj = 1, jpj 
    466464            DO ji = 1, jpi 
    467                IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    468                   ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    469                   old_v_i(ji,jj,jl)     = d_v_i_trp(ji,jj,jl) 
    470                   d_v_i_trp(ji,jj,jl)   = 0._wp 
    471                   old_a_i(ji,jj,jl)     = d_a_i_trp(ji,jj,jl) 
    472                   d_a_i_trp(ji,jj,jl)   = 0._wp 
    473                   old_v_s(ji,jj,jl)     = d_v_s_trp(ji,jj,jl) 
    474                   d_v_s_trp(ji,jj,jl)   = 0._wp 
    475                   old_e_s(ji,jj,1,jl)  = d_e_s_trp(ji,jj,1,jl) 
    476                   d_e_s_trp(ji,jj,1,jl) = 0._wp 
    477                   old_oa_i(ji,jj,jl)    = d_oa_i_trp(ji,jj,jl) 
    478                   d_oa_i_trp(ji,jj,jl)  = 0._wp 
    479                   IF(  num_sal == 2  .OR.  num_sal == 4  )   old_smv_i(ji,jj,jl)  = d_smv_i_trp(ji,jj,jl) 
    480                   d_smv_i_trp(ji,jj,jl) = 0._wp 
     465               IF(  old_v_i  (ji,jj,jl) < epsi06 .AND. & 
     466                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
     467                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
     468                  d_v_i_trp (ji,jj,jl)   = 0._wp 
     469                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
     470                  d_a_i_trp (ji,jj,jl)   = 0._wp 
     471                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
     472                  d_v_s_trp (ji,jj,jl)   = 0._wp 
     473                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
     474                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
     475                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
     476                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
     477                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
     478                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
    481479               ENDIF 
    482480            END DO 
    483481         END DO 
    484482      END DO 
    485  
     483      ! 
    486484      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    487485      ! 
    … …  
    612610                  ! present 
    613611                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    614                      + strength(ji-1,jj) * tms(ji-1,jj) &   
    615                      + strength(ji+1,jj) * tms(ji+1,jj) &   
    616                      + strength(ji,jj-1) * tms(ji,jj-1) &   
    617                      + strength(ji,jj+1) * tms(ji,jj+1)     
     612                     &          + strength(ji-1,jj) * tms(ji-1,jj) &   
     613                     &          + strength(ji+1,jj) * tms(ji+1,jj) &   
     614                     &          + strength(ji,jj-1) * tms(ji,jj-1) &   
     615                     &          + strength(ji,jj+1) * tms(ji,jj+1)     
    618616 
    619617                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 
    620618                  zworka(ji,jj) = zworka(ji,jj) / zw1 
    621619               ELSE 
    622                   zworka(ji,jj) = 0.0 
     620                  zworka(ji,jj) = 0._wp 
    623621               ENDIF 
    624622            END DO 
    … …  
    10481046         DO jj = 1, jpj 
    10491047            DO ji = 1, jpi 
    1050                IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0       & 
    1051                   .AND. closing_gross(ji,jj) > 0.0) THEN 
     1048               IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       & 
     1049                  .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    10521050                  icells = icells + 1 
    10531051                  indxi(icells) = ji 
    … …  
    11301128            ! Salinity 
    11311129            !------------- 
    1132             smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of water frozen in voids 
     1130            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids 
    11331131 
    11341132            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    … …  
    11371135             
    11381136            !                                                             ! excess of salt is flushed into the ocean 
    1139             fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 
    1140  
     1137            fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
     1138 
     1139            rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0   ! increase in ice volume du to seawater frozen in voids 
     1140             
    11411141            !------------------------------------             
    11421142            ! 3.6 Increment ridging diagnostics 
    … …  
    11481148            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    11491149            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1150             diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice 
    1151             opening    (ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice 
     1150            diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
     1151            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    11521152 
    11531153            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 
    … …  
    11561156            ! 3.7 Put the snow somewhere in the ocean 
    11571157            !------------------------------------------             
    1158  
    11591158            !  Place part of the snow lost by ridging into the ocean.  
    11601159            !  Note that esnow_mlt < 0; the ocean must cool to melt snow. 
    … …  
    11791178            !           ij looping 1-icells 
    11801179 
    1181             dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
     1180            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
    11821181            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 
    1183  
    11841182 
    11851183         END DO                 ! ij 
    … …  
    12111209 
    12121210               ! heat flux 
    1213                fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / rdt_ice 
     1211               fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 
    12141212 
    12151213               ! Correct dimensions to avoid big values 
    … …  
    12751273               ! Transfer area, volume, and energy accordingly. 
    12761274 
    1277                IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR.        & 
    1278                   hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN 
    1279                   hL = 0.0 
    1280                   hR = 0.0 
     1275               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR.        & 
     1276                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
     1277                  hL = 0._wp 
     1278                  hR = 0._wp 
    12811279               ELSE 
    1282                   hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1)) 
    1283                   hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2)) 
     1280                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 
     1281                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   ) 
    12841282               ENDIF 
    12851283 
    12861284               ! fraction of ridged ice area and volume going to n2 
    1287                farea = (hR-hL) / dhr(ji,jj)  
    1288                fvol(ji,jj) = (hR*hR - hL*hL) / dhr2(ji,jj) 
    1289  
    1290                a_i  (ji,jj,jl2)   = a_i  (ji,jj,jl2)  + ardg2 (ji,jj) * farea 
    1291                v_i  (ji,jj,jl2)   = v_i  (ji,jj,jl2)  + vrdg2 (ji,jj) * fvol(ji,jj) 
    1292                v_s  (ji,jj,jl2)   = v_s  (ji,jj,jl2)  + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
     1285               farea = ( hR - hL ) / dhr(ji,jj)  
     1286               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj) 
     1287 
     1288               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea 
     1289               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 
     1290               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
    12931291               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
    1294                smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2)  + srdg2 (ji,jj) * fvol(ji,jj) 
    1295                oa_i (ji,jj,jl2)   = oa_i (ji,jj,jl2)  + oirdg2(ji,jj) * farea 
     1292               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 
     1293               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    12961294 
    12971295            END DO ! ij 
    … …  
    13171315               ! Compute the fraction of rafted ice area and volume going to  
    13181316               ! thickness category jl2, transfer area, volume, and energy accordingly. 
    1319  
    1320                IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1321                   hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
    1322                   a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) 
    1323                   v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) 
    1324                   v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + vsrft(ji,jj)*fsnowrft 
    1325                   e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + esrft(ji,jj)*fsnowrft 
    1326                   smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + smrft(ji,jj)     
    1327                   oa_i(ji,jj,jl2)  = oa_i(ji,jj,jl2) + oirft2(ji,jj)     
     1317               ! 
     1318               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND.        & 
     1319                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     1320                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj) 
     1321                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj) 
     1322                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft 
     1323                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft 
     1324                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
     1325                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
    13281326               ENDIF ! hraft 
    1329  
     1327               ! 
    13301328            END DO ! ij 
    13311329 
    … …  
    13361334                  ji = indxi(ij) 
    13371335                  jj = indxj(ij) 
    1338                   IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1339                      hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
     1336                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)  .AND.        & 
     1337                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN 
    13401338                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    13411339                  ENDIF 
    … …  
    15041502            DO jj = 1 , jpj 
    15051503               DO ji = 1 , jpi 
    1506 !!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 
     1504!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice 
    15071505!!gm                  xtmp = xtmp * unit_fac 
    15081506                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    … …  
    15241522               ! fluxes are positive to the ocean 
    15251523               ! here the flux has to be negative for the ocean 
    1526 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 
     1524!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice 
    15271525               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    15281526 
    1529 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB   ??????? 
     1527!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    15301528 
    15311529               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
    … …  
    15401538 
    15411539               !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj)                  ) * &  
    1542                !                               rhosn * v_s(ji,jj,jl) / rdt_ice 
     1540               !                               rhosn * v_s(ji,jj,jl) * r1_rdtice 
    15431541 
    15441542               !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
    1545                !                               rhoic * v_i(ji,jj,jl) / rdt_ice 
     1543               !                               rhoic * v_i(ji,jj,jl) * r1_rdtice 
    15461544 
    15471545               !           sfx (i,j)      = sfx (i,j)      + xtmp 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r3294 r3517  
    101101 
    102102      !- Trend terms 
    103       d_a_i_thd (:,:,:)  = a_i(:,:,:)   - old_a_i(:,:,:)  
    104       d_v_s_thd (:,:,:)  = v_s(:,:,:)   - old_v_s(:,:,:) 
    105       d_v_i_thd (:,:,:)  = v_i(:,:,:)   - old_v_i(:,:,:)   
     103      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
     104      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
     105      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
    106106      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
    107107      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    108108 
    109109      d_smv_i_thd(:,:,:) = 0._wp 
    110       IF( num_sal == 2 .OR. num_sal == 4 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     110      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    111111 
    112112      IF(ln_ctl) THEN   ! Control print 
    … …  
    143143 
    144144      !- Recover Old values 
    145       a_i(:,:,:)   = old_a_i (:,:,:) 
    146       v_s(:,:,:)   = old_v_s (:,:,:) 
    147       v_i(:,:,:)   = old_v_i (:,:,:) 
    148       e_s(:,:,:,:) = old_e_s (:,:,:,:) 
    149       e_i(:,:,:,:) = old_e_i (:,:,:,:) 
    150       ! 
    151       IF( num_sal == 2 .OR. num_sal == 4 )   smv_i(:,:,:)       = old_smv_i (:,:,:) 
     145      a_i(:,:,:)   = old_a_i(:,:,:) 
     146      v_s(:,:,:)   = old_v_s(:,:,:) 
     147      v_i(:,:,:)   = old_v_i(:,:,:) 
     148      e_s(:,:,:,:) = old_e_s(:,:,:,:) 
     149      e_i(:,:,:,:) = old_e_i(:,:,:,:) 
     150      ! 
     151      IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    152152      ! 
    153153   END SUBROUTINE lim_itd_th 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r3294 r3517  
    402402                     zsmax = 4.5_wp 
    403403                     zsmin = 3.5_wp 
    404                      IF( sm_i(ji,jj,jl) .LT. zsmin ) THEN 
     404                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN 
    405405                        zalpha = 1._wp 
    406                      ELSEIF( sm_i(ji,jj,jl) .LT.zsmax ) THEN 
     406                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN 
    407407                        zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin ) 
    408408                     ELSE 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r3488 r3517  
    99   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea 
    1010   !!                 !                  + simplification of the ice-ocean stress calculation 
    11    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     12   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim3 
    … …  
    4445   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    4546 
    46    REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice  
    4747   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values 
    4848   REAL(wp)  ::   rzero  = 0._wp     
    … …  
    8888      !!              - qns     : sea heat flux: non solar 
    8989      !!              - emp     : freshwater budget: volume flux  
    90       !!              - sfx     : freshwater budget: concentration/dillution  
     90      !!              - sfx     : salt flux  
    9191      !!              - fr_i    : ice fraction 
    9292      !!              - tn_ice  : sea-ice surface temperature 
    … …  
    102102      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    103103      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    104       REAL(wp) ::   zinda, zfons, zpme              ! local scalars 
    105       REAL(wp), POINTER, DIMENSION(:,:) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
    106       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace 
     104      REAL(wp) ::   zinda, zemp, zemp_snow, zfmm    ! local scalars 
     105      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfcm1 , zfcm2   ! solar/non solar heat fluxes 
     106      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
    107107      !!--------------------------------------------------------------------- 
    108108       
    … …  
    153153               !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used 
    154154               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
    155                &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                  * r1_rdtice   & 
    156                &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!! 
     155               &           + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   & 
     156               &           + fhmec(ji,jj)     & ! new contribution due to snow melt when ridging!! 
    157157               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation 
    158158               &           + fheat_res(ji,jj) 
    159             ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean 
    160             !         computed in limthd_zdf.F90 
     159            ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean computed in limthd_zdf.F90 
    161160            ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t 
    162161            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok) 
    … …  
    167166            ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead 
    168167 
    169             IF ( num_sal == 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + & 
    170                fhbri(ji,jj) ! new contribution due to brine drainage  
    171  
    172             ! bottom radiative component is sent to the computation of the 
    173             ! oceanic heat flux 
     168            IF( num_sal == 2 )   zfcm2(ji,jj) = zfcm2(ji,jj) + fhbri(ji,jj)    ! add contribution due to brine drainage  
     169 
     170            ! bottom radiative component is sent to the computation of the oceanic heat flux 
    174171            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)      
    175172 
    … …  
    179176            !                           ! fdtcn : turbulent oceanic heat flux 
    180177 
    181             !!gm   this IF prevents the vertorisation of the whole loop 
     178!!gm   this IF prevents the vertorisation of the whole loop 
    182179            IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    183180               WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    … …  
    208205               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    209206            ENDIF 
    210             !!gm   end 
     207!!gm   end 
    211208         END DO 
    212209      END DO 
    … …  
    229226 
    230227            !  computing freshwater exchanges at the ice/ocean interface 
    231             zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
    232                &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
    233                &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
    234                &   - rdm_snw(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting  
    235                &   + fmmec(ji,jj)                                         ! snow falling when ridging 
    236  
    237  
    238             !  computing salt exchanges at the ice/ocean interface 
    239             !  sice should be the same as computed with the ice model 
    240             zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdm_ice(ji,jj) * r1_rdtice  
    241             ! SOCE 
    242             zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdm_ice(ji,jj) * r1_rdtice 
    243  
    244             !CT useless            !  salt flux for constant salinity 
    245             !CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj) 
    246             !  salt flux for variable salinity 
    247             zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
    248             !  correcting brine and salt fluxes 
    249             fsbri(ji,jj)      =  zinda*fsbri(ji,jj) 
    250             !  converting the salt fluxes from ice to a freshwater flux from ocean 
    251             fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 
    252             fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 ) 
    253             fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 ) 
    254             fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 
    255  
    256             !  freshwater mass exchange (positive to the ice, negative for the ocean ?) 
    257             !  actually it's a salt flux (so it's minus freshwater flux) 
    258             !  if sea ice grows, zfons is positive, fsalt also 
    259             !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN 
    260             !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1] 
    261  
    262             emp(ji,jj) = - zpme  
     228            zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
     229               &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
     230               &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
     231               &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     232 
     233            ! mass flux at the ocean/ice interface (sea ice fraction) 
     234            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean 
     235            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting   
     236 
     237            emp(ji,jj) = zemp + zemp_snw + zfmm  ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     238             
     239            !  correcting brine salt fluxes   (zinda = 1  if pfrld=1 , =0 otherwise) 
     240            zinda        = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
     241            fsbri(ji,jj) = zinda * fsbri(ji,jj) 
    263242         END DO 
    264243      END DO 
    265244 
     245      !------------------------------------------! 
     246      !      salt flux at the ocean surface      ! 
     247      !------------------------------------------! 
     248 
    266249      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux 
    267          sfx (:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
     250         sfx (:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + fsbri(:,:) 
    268251      ELSE                         ! constant ice salinity: 
    269          sfx (:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
     252         sfx (:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) 
    270253      ENDIF 
    271254      !-----------------------------------------------! 
    … …  
    277260         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )  
    278261         !                                                      ! time evolution of snow+ice mass 
    279          snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice 
     262         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice 
    280263      ENDIF 
    281264 
    … …  
    403386      !                                      ! allocate lim_sbc array 
    404387      IF( lim_sbc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' ) 
    405       ! 
    406       r1_rdtice = 1. / rdt_ice 
    407388      ! 
    408389      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r3419 r3517  
    110110                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i 
    111111                  !0 if no ice and 1 if yes 
    112                   zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) )  
     112                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) )  
    113113                  !convert units ! very important that this line is here 
    114114                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    122122                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s 
    123123                  !0 if no ice and 1 if yes 
    124                   zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) )  
     124                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) )  
    125125                  !convert units ! very important that this line is here 
    126126                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    284284            CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb), fbif       , jpi, jpj, npb(1:nbpb) ) 
    285285            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif      , jpi, jpj, npb(1:nbpb) ) 
    286             CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb), rdm_ice    , jpi, jpj, npb(1:nbpb) ) 
    287             CALL tab_2d_1d( nbpb, rdmsnif_1d (1:nbpb), rdm_snw    , jpi, jpj, npb(1:nbpb) ) 
     286            CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice    , jpi, jpj, npb(1:nbpb) ) 
     287            CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw    , jpi, jpj, npb(1:nbpb) ) 
    288288            CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi      , jpi, jpj, npb(1:nbpb) ) 
    289289            CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq    , jpi, jpj, npb(1:nbpb) ) 
    … …  
    352352            CALL tab_1d_2d( nbpb, qldif  , npb, qldif_1d  (1:nbpb), jpi, jpj ) 
    353353            CALL tab_1d_2d( nbpb, qfvbq  , npb, qfvbq_1d  (1:nbpb), jpi, jpj ) 
    354             CALL tab_1d_2d( nbpb, rdm_ice, npb, rdmicif_1d(1:nbpb), jpi, jpj ) 
    355             CALL tab_1d_2d( nbpb, rdm_snw, npb, rdmsnif_1d(1:nbpb), jpi, jpj ) 
     354            CALL tab_1d_2d( nbpb, rdm_ice, npb, rdm_ice_1d(1:nbpb), jpi, jpj ) 
     355            CALL tab_1d_2d( nbpb, rdm_snw, npb, rdm_snw_1d(1:nbpb), jpi, jpj ) 
    356356            CALL tab_1d_2d( nbpb, dmgwi  , npb, dmgwi_1d  (1:nbpb), jpi, jpj ) 
    357357            CALL tab_1d_2d( nbpb, rdvosif, npb, dvsbq_1d  (1:nbpb), jpi, jpj ) 
    … …  
    419419      !-------------------------------------------- 
    420420      d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes  
    421       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0 
     421      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    422422 
    423423      IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
    … …  
    488488      ! 
    489489      IF(lwp) WRITE(numout,*) ' lim_thd_glohec ' 
    490       IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice 
    491       IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice 
    492       IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) / rdt_ice 
     490      IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
     491      IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
     492      IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
    493493      ! 
    494494   END SUBROUTINE lim_thd_glohec 
    … …  
    538538      !-------------------- 
    539539      DO ji = kideb, kiut 
    540          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     540         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    541541      END DO 
    542542 
    … …  
    597597            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 
    598598            WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 
    599             WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) / rdt_ice 
     599            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) * r1_rdtice 
    600600            WRITE(numout,*) ' Fdt        : ', sum_fluxq(ji,jl) 
    601601            WRITE(numout,*) 
    … …  
    631631            WRITE(numout,*) 
    632632            WRITE(numout,*) ' Layer by layer ... ' 
    633             WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 
     633            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 
    634634            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) - fc_s(ji,0) 
    635635            DO jk = 1, nlay_i 
    636636               WRITE(numout,*) ' layer  : ', jk 
    637                WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice   
     637               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) * r1_rdtice   
    638638               WRITE(numout,*) ' radab  : ', radab(ji,jk) 
    639639               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) - fc_i(ji,jk-1) 
    … …  
    681681         fatm      (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji)                       ! total heat flux 
    682682         sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(zji,zjj,jl)  
    683          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     683         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    684684      END DO 
    685685 
    … …  
    688688      !-------------------- 
    689689      DO ji = kideb, kiut 
    690          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     690         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    691691      END DO 
    692692 
    … …  
    722722            WRITE(numout,*) ' * ' 
    723723            WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl) 
    724             WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) / rdt_ice 
    725             WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice 
    726             WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 
     724            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) * r1_rdtice 
     725            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) * r1_rdtice 
     726            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 
    727727            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 
    728728            WRITE(numout,*) ' * ' 
    … …  
    734734            WRITE(numout,*) ' * ' 
    735735            WRITE(numout,*) ' Heat contents --- : ' 
    736             WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice 
    737             WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice 
    738             WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice 
    739             WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice 
    740             WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice 
    741             WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice 
     736            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) * r1_rdtice 
     737            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) * r1_rdtice 
     738            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) * r1_rdtice 
     739            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) * r1_rdtice 
     740            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) * r1_rdtice 
     741            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) * r1_rdtice 
    742742            WRITE(numout,*) ' * ' 
    743743            WRITE(numout,*) ' Ice variables --- : ' 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r3294 r3517  
    77   !!                 ! 2005-06 (M. Vancoppenolle) 3D version  
    88   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 
    9    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     9   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     10   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes  
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim3 
    … …  
    7172      !!  
    7273      INTEGER  ::   ji , jk        ! dummy loop indices 
    73       INTEGER  ::   zji, zjj       ! 2D corresponding indices to ji 
     74      INTEGER  ::   ii, ij       ! 2D corresponding indices to ji 
    7475      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7576      INTEGER  ::   isnowic        ! snow ice formation not 
    … …  
    145146      ! 
    146147      DO ji = kideb, kiut 
    147          isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 
    148          ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt 
    149          z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
    150          z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 
     148         isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  ) ) 
     149         ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 
     150         z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
     151         z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 
    151152         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 
    152153      END DO ! ji 
    … …  
    240241         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji) 
    241242         ! If snow is still present zhn = 1, else zhn = 0 
    242          zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 
     243         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew ) ) 
    243244         ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
    244245         ! Volume and mass variations of snow 
    245          dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 
     246         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 
    246247         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) ) 
    247          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji) 
     248         rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
    248249      END DO ! ji 
    249250 
    … …  
    253254      DO ji = kideb, kiut  
    254255         dh_i_surf(ji) =  0._wp 
    255          z_f_surf (ji) =  zqfont_su(ji) / rdt_ice ! heat conservation test 
     256         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice  ! heat conservation test 
    256257         zdq_i    (ji) =  0._wp 
    257258      END DO ! ji 
    … …  
    262263            zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 
    263264            !                                                    ! recompute heat available 
    264             zqfont_su(ji)    = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
     265            zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
    265266            !                                                    ! melt of layer jk cannot be higher than its thickness 
    266267            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 
    267268            !                                                    ! update surface melt 
    268             dh_i_surf(ji)    = dh_i_surf(ji) + zdeltah(ji,jk)  
     269            dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk)  
    269270            !                                                    ! for energy conservation 
    270             zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
     271            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    271272            ! 
    272             ! contribution to ice-ocean salt flux  
    273             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    274             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    275             zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji)    & 
    276                &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice  
     273            !                                                    ! contribution to ice-ocean salt flux  
     274            zfsalt_melt(ji)  = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
    277275         END DO 
    278276      END DO 
    … …  
    290288            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN! 
    291289               WRITE(numout,*) ' ALERTE heat loss for surface melt ' 
    292                WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 
     290               WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl 
    293291               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji) 
    294292               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji) 
    … …  
    299297               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji) 
    300298               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji) 
    301                WRITE(numout,*) ' sss_m        : ', sss_m(zji,zjj) 
     299               WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij) 
    302300            ENDIF 
    303301         END DO 
    … …  
    338336         DO ji = kideb, kiut 
    339337            ! In case of disparition of the snow, we have to update the snow temperatures 
    340             zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 
     338            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) ) ) 
    341339            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt 
    342340            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) 
    … …  
    358356      ! 4.1 Basal growth - (a) salinity not varying in time  
    359357      !----------------------------------------------------- 
    360       IF(  num_sal /= 2  .AND.  num_sal /= 4  ) THEN 
    361          DO ji = kideb, kiut 
    362             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0  ) THEN 
     358      IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time 
     359         DO ji = kideb, kiut 
     360            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN 
    363361               s_i_new(ji)         =  sm_i_b(ji) 
    364362               ! Melting point in K 
    … …  
    371369                  &                           - rcp  * ( ztmelts - rtt )                                 ) 
    372370               ! Basal growth rate = - F*dt / q 
    373                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     371               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
    374372            ENDIF 
    375373         END DO 
    … …  
    379377      ! 4.1 Basal growth - (b) salinity varying in time  
    380378      !------------------------------------------------- 
    381       IF(  num_sal == 2 .OR.  num_sal == 4  ) THEN 
    382          ! the growth rate (dh_i_bott) is function of the new ice 
    383          ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice  
    384          ! salinity (snewice). snewice depends on dh_i_bott 
    385          ! it converges quickly, so, no problem 
     379      IF(  num_sal == 2  ) THEN 
     380         ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).  
     381         ! q_i_b depends on the new ice salinity (snewice).  
     382         ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 
    386383         ! See Vancoppenolle et al., OM08 for more info on this 
    387384 
    … …  
    394391            DO ji = kideb, kiut 
    395392               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
    396                   zji = MOD( npb(ji) - 1, jpi ) + 1 
    397                   zjj = ( npb(ji) - 1 ) / jpi + 1 
     393                  ii = MOD( npb(ji) - 1, jpi ) + 1 
     394                  ij = ( npb(ji) - 1 ) / jpi + 1 
    398395                  ! Melting point in K 
    399396                  ztmelts             =   - tmut * s_i_new(ji) + rtt  
    … …  
    408405                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    409406                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    410                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) ) 
     407                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
    411408                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )  
    412409                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    … …  
    414411                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    415412                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    416                   zds         = zfracs * sss_m(zji,zjj) - s_i_new(ji) 
    417                   s_i_new(ji) = zfracs * sss_m(zji,zjj) 
     413                  zds         = zfracs * sss_m(ii,ij) - s_i_new(ji) 
     414                  s_i_new(ji) = zfracs * sss_m(ii,ij) 
    418415               ENDIF ! fc_bo_i 
    419416            END DO ! ji 
    … …  
    432429                  &                            - rcp * ( ztmelts - rtt )                                    ) 
    433430               ! Basal growth rate = - F*dt / q 
    434                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     431               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    435432               ! Salinity update 
    436433               ! entrapment during bottom growth 
    … …  
    453450            s_i_new(ji)   =  s_i_b(ji,nlay_i) 
    454451            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 
    455             zfbase(ji)    =  zqfont_bo(ji) / rdt_ice     ! heat conservation test 
     452            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test 
    456453            zdq_i(ji)     =  0._wp 
    457454            dh_i_bott(ji) =  0._wp 
    … …  
    461458      DO jk = nlay_i, 1, -1 
    462459         DO ji = kideb, kiut 
    463             IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN 
    464                ztmelts            =   - tmut * s_i_b(ji,jk) + rtt  
    465                IF( t_i_b(ji,jk) >= ztmelts ) THEN 
    466                   zdeltah(ji,jk)  = - zh_i(ji) 
    467                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk) 
    468                   zinnermelt(ji)   = 1._wp 
    469                ELSE  ! normal ablation 
    470                   zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk) 
    471                   zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
    472                   zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
    473                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk) 
    474                   zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
    475                   ! contribution to salt flux 
    476                   zji             = MOD( npb(ji) - 1, jpi ) + 1 
    477                   zjj             = ( npb(ji) - 1 ) / jpi + 1 
    478                   zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji)   ) * a_i_b(ji)   & 
    479                      &                              * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice  
     460            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN 
     461               ztmelts = - tmut * s_i_b(ji,jk) + rtt  
     462               IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed 
     463                  zdeltah   (ji,jk) = - zh_i(ji) 
     464                  dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
     465                  zinnermelt(ji   ) = 1._wp 
     466               ELSE                                  ! normal ablation 
     467                  zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 
     468                  zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
     469                  zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
     470                  dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
     471                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    480472               ENDIF 
     473               ! contribution to salt flux 
     474               zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
    481475            ENDIF 
    482476         END DO ! ji 
    … …  
    493487               ENDIF 
    494488               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN 
    495                   WRITE(numout,*) ' ALERTE heat loss for basal melt : zji, zjj, jl :', zji, zjj, jl 
     489                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl 
    496490                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji) 
    497491                  WRITE(numout,*) ' zfbase    : ', zfbase(ji) 
    … …  
    502496                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji) 
    503497                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji) 
    504                   WRITE(numout,*) ' sss_m     : ', sss_m(zji,zjj) 
     498                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij) 
    505499                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    506500                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) 
    … …  
    531525         !                     ! excessive energy is sent to lateral ablation 
    532526         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   & 
    533             &                          * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 
     527            &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 
    534528         dh_i_bott(ji)  = zdhbf 
    535529         !                     !since ice volume is only used for outputs, we keep it global for all categories 
    … …  
    538532         zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    539533         !                     ! diagnostic ( bottom ice growth ) 
    540          zji = MOD( npb(ji) - 1, jpi ) + 1 
    541          zjj = ( npb(ji) - 1 ) / jpi + 1 
    542          diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 
    543          diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice 
    544          diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 
     534         ii = MOD( npb(ji) - 1, jpi ) + 1 
     535         ij = ( npb(ji) - 1 ) / jpi + 1 
     536         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     537         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
     538         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    545539      END DO 
    546540 
    … …  
    548542      ! 5.2 More than available ice melts 
    549543      !----------------------------------- 
    550       ! then heat applied minus heat content at previous time step 
    551       ! should equal heat remaining  
     544      ! then heat applied minus heat content at previous time step should equal heat remaining  
    552545      ! 
    553546      DO ji = kideb, kiut 
    554547         ! Adapt the remaining energy if too much ice melts 
    555548         !-------------------------------------------------- 
    556          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
     549         zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    557550         ! 0 if no more ice 
    558551         zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0 
    … …  
    562555         ! If snow remains, energy is used to melt snow 
    563556         zhni =  ht_s_b(ji)      ! snow depth at previous time step 
    564          zihg =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow  
     557         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow  
    565558 
    566559         ! energy of melting of remaining snow 
    567560         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 
    568561         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 
    569          zhnfi          =  zhni + zdhnm 
     562         zhnfi     =  zhni + zdhnm 
    570563         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 
    571564         ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    … …  
    581574         ! 
    582575         !                                              ! mass variation cumulated over category 
    583          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s                     ! snow  
    584          rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i                     ! ice  
     576         rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
     577         rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
    585578 
    586579         ! Remaining heat to the ocean  
    587580         !--------------------------------- 
    588          focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt 
    589  
    590       END DO 
    591  
    592       ftotal_fin (:) = zfdt_final(:)  / rdt_ice 
     581         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt 
     582 
     583      END DO 
     584 
     585      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice 
    593586 
    594587      !--------------------------- 
    … …  
    596589      !--------------------------- 
    597590      DO ji = kideb, kiut 
    598          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   !1 if ice 
    599  
     591         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
     592         ! 
    600593         ! Salt flux 
    601          zji = MOD( npb(ji) - 1, jpi ) + 1 
    602          zjj = ( npb(ji) - 1 ) / jpi + 1 
    603          ! new lines 
    604          IF( num_sal == 4 ) THEN 
    605             fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                & 
    606                &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal   ) / rdt_ice 
    607          ELSE 
    608             fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                & 
    609                &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 
    610          ENDIF 
     594         fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                & 
     595            &                        - (1.0 - zihgnew) * zfmass_i   (ji) * sm_i_b(ji) ) * r1_rdtice 
     596         ! 
    611597         ! Heat flux 
    612598         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 
    613          ! excessive total ablation energy (focea) sent to the ocean 
     599         ! excessive total  ablation energy (focea) sent to the ocean 
    614600         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 
    615601 
    616          zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 
    617          ! equals 0 if ht_i = 0, 1 if ht_i gt 0 
     602         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0 
    618603         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji) 
    619          qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji)    * a_i_b(ji) * rdt_ice   & 
     604         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   & 
    620605            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice 
    621606      END DO  ! ji 
    … …  
    656641         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 
    657642 
    658          rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
    659          rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
     643         rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
     644         rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
    660645 
    661646         !        Equivalent salt flux (1) Snow-ice formation component 
    662647         !        ----------------------------------------------------- 
    663          zji = MOD( npb(ji) - 1, jpi ) + 1 
    664          zjj =    ( npb(ji) - 1 ) / jpi + 1 
    665  
    666          IF( num_sal /= 2 ) THEN   ;   zsm_snowice = sm_i_b(ji) 
    667          ELSE                      ;   zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj)  
     648         ii = MOD( npb(ji) - 1, jpi ) + 1 
     649         ij =    ( npb(ji) - 1 ) / jpi + 1 
     650 
     651         IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 
     652         ELSE                      ;   zsm_snowice = sm_i_b(ji)   
    668653         ENDIF 
    669          IF( num_sal == 4 ) THEN 
    670             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    ) * a_i_b(ji)   & 
    671                &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 
    672          ELSE 
    673             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji)   & 
    674                &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 
    675          ENDIF 
     654         fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice  * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     655         ! 
    676656         ! entrapment during snow ice formation 
    677657         i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 
    678658         isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch 
    679          IF(  num_sal == 2  .OR.  num_sal == 4  )   & 
    680             dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 
    681             &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   & 
    682             &               - sm_i_b(ji) ) * isnowic      
     659         IF(  num_sal == 2  )   & 
     660            dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    & 
     661            &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   & 
     662            &                  - sm_i_b(ji) ) * isnowic      
    683663 
    684664         !  Actualize new snow and ice thickness. 
    … …  
    690670 
    691671         ! diagnostic ( snow ice growth ) 
    692          zji = MOD( npb(ji) - 1, jpi ) + 1 
    693          zjj =    ( npb(ji) - 1 ) / jpi + 1 
    694          diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice 
     672         ii = MOD( npb(ji) - 1, jpi ) + 1 
     673         ij =    ( npb(ji) - 1 ) / jpi + 1 
     674         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    695675         ! 
    696676      END DO !ji 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r3351 r3517  
    147147      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis 
    148148      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms 
    149       !!------------------------------------------------------------------ 
    150        
     149      !!------------------------------------------------------------------      
    151150      !  
    152151      !------------------------------------------------------------------------------! 
    … …  
    156155      DO ji = kideb , kiut 
    157156         ! is there snow or not 
    158          isnow(ji)= INT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
     157         isnow(ji)= INT(  1._wp - MAX(  0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
    159158         ! surface temperature of fusion 
    160159!!gm ???  ztfs(ji) = rtt !!!???? 
    … …  
    201200      DO ji = kideb , kiut 
    202201         ! switches 
    203          isnow(ji) = INT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
     202         isnow(ji) = INT(  1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
    204203         ! hs > 0, isnow = 1 
    205204         zhsu (ji) = hnzst  ! threshold for the computation of i0 
    … …  
    262261      ! just to check energy conservation 
    263262      DO ji = kideb, kiut 
    264          ii                = MOD( npb(ji) - 1, jpi ) + 1 
    265          ij                = ( npb(ji) - 1 ) / jpi + 1 
     263         ii = MOD( npb(ji) - 1 , jpi ) + 1 
     264         ij =    ( npb(ji) - 1 ) / jpi + 1 
    266265         fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i) 
    267266      END DO 
    … …  
    273272         END DO 
    274273      END DO 
    275  
    276274 
    277275      ! 
    … …  
    662660 
    663661            ! surface temperature 
    664             isnow(ji)     = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) 
     662            isnow(ji)     = INT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  ) 
    665663            ztsuoldit(ji) = t_su_b(ji) 
    666             IF (t_su_b(ji) .LT. ztfs(ji)) & 
     664            IF( t_su_b(ji) < ztfs(ji) )  & 
    667665               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1)   & 
    668666               &          + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r3294 r3517  
    408408      IF ( con_i ) THEN 
    409409         DO ji = kideb, kiut 
    410             IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN 
     410            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  > 1.0e-6 ) THEN 
    411411               zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    412412               zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    413413               WRITE(numout,*) ' violation of heat conservation : ',             & 
    414                   ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice 
     414                  ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 
    415415               WRITE(numout,*) ' ji, jj   : ', zji, zjj 
    416416               WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji) 
    417                WRITE(numout,*) ' zqts_in  : ', zqts_in(ji) / rdt_ice 
    418                WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice 
     417               WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice 
     418               WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice 
    419419               WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) 
    420420               WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 
    … …  
    526526         ! bottom formation temperature 
    527527         ztform = t_i_b(ji,nlay_i) 
    528          IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) 
     528         IF(  num_sal == 2  )  ztform = t_bo_b(ji) 
    529529         qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice 
    530530            &              + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice 
    … …  
    622622      ! 
    623623      DO ji = kideb, kiut 
    624          IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN 
     624         IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  > 1.0e-6 ) THEN 
    625625            zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    626626            zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    627             WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice 
     627            WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 
    628628            WRITE(numout,*) ' ji, jj   : ', zji, zjj 
    629629            WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji) 
    630             WRITE(numout,*) ' zqti_in  : ', zqti_in(ji) / rdt_ice 
    631             WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) / rdt_ice 
     630            WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice 
     631            WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 
    632632            WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 
    633633            WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3294 r3517  
    143143      ! 1) Conservation check and changes in each ice category 
    144144      !------------------------------------------------------------------------------! 
    145       IF ( con_i ) THEN 
    146          CALL lim_column_sum (jpl, v_i, vt_i_init) 
    147          CALL lim_column_sum (jpl, v_s, vt_s_init) 
    148          CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init) 
    149          CALL lim_column_sum (jpl,  e_s(:,:,1,:) , et_s_init) 
     145      IF( con_i ) THEN 
     146         CALL lim_column_sum        ( jpl, v_i          , vt_i_init) 
     147         CALL lim_column_sum        ( jpl, v_s          , vt_s_init) 
     148         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 
     149         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init) 
    150150      ENDIF 
    151151 
    … …  
    158158               DO ji = 1, jpi 
    159159                  !Energy of melting q(S,T) [J.m-3] 
    160                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / & 
    161                      MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * & 
    162                      nlay_i 
    163                   zindb      = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes 
    164                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb 
     160                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * nlay_i 
     161                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
     162                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
    165163               END DO 
    166164            END DO 
    … …  
    182180      !  
    183181 
    184       zvrel(:,:) = 0.0 
     182      zvrel(:,:) = 0._wp 
    185183 
    186184      ! Default new ice thickness  
    187       DO jj = 1, jpj 
    188          DO ji = 1, jpi 
    189             hicol(ji,jj) = hiccrit(1) 
    190          END DO 
    191       END DO 
    192  
    193       IF (fraz_swi.eq.1.0) THEN 
     185      hicol(:,:) = hiccrit(1) 
     186 
     187      IF( fraz_swi == 1._wp ) THEN 
    194188 
    195189         !-------------------- 
    196190         ! Physical constants 
    197191         !-------------------- 
    198          hicol(:,:) = 0.0 
     192         hicol(:,:) = 0._wp 
    199193 
    200194         zhicrit = 0.04 ! frazil ice thickness 
    … …  
    306300      IF ( nbpac > 0 ) THEN 
    307301 
    308          CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         ,       & 
    309             jpi, jpj, npac(1:nbpac) ) 
     302         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) ) 
    310303         DO jl = 1, jpl 
    311             CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  ,       & 
    312                jpi, jpj, npac(1:nbpac) ) 
    313             CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  ,       & 
    314                jpi, jpj, npac(1:nbpac) ) 
    315             CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) ,       & 
    316                jpi, jpj, npac(1:nbpac) ) 
    317             CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl),       & 
    318                jpi, jpj, npac(1:nbpac) ) 
     304            CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) ) 
     305            CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) ) 
     306            CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , jpi, jpj, npac(1:nbpac) ) 
     307            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    319308            DO jk = 1, nlay_i 
    320                CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , & 
    321                   jpi, jpj, npac(1:nbpac) ) 
     309               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    322310            END DO ! jk 
    323311         END DO ! jl 
    324312 
    325          CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif ,              & 
    326             jpi, jpj, npac(1:nbpac) ) 
    327          CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif ,              & 
    328             jpi, jpj, npac(1:nbpac) ) 
    329          CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  ,              & 
    330             jpi, jpj, npac(1:nbpac) ) 
    331          CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv ,              & 
    332             jpi, jpj, npac(1:nbpac) ) 
    333          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol ,              & 
    334             jpi, jpj, npac(1:nbpac) ) 
    335          CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel ,              & 
    336             jpi, jpj, npac(1:nbpac) ) 
     313         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif , jpi, jpj, npac(1:nbpac) ) 
     314         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif , jpi, jpj, npac(1:nbpac) ) 
     315         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  , jpi, jpj, npac(1:nbpac) ) 
     316         CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv , jpi, jpj, npac(1:nbpac) ) 
     317         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol , jpi, jpj, npac(1:nbpac) ) 
     318         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel , jpi, jpj, npac(1:nbpac) ) 
    337319 
    338320         !------------------------------------------------------------------------------! 
    … …  
    344326         !---------------------- 
    345327         DO ji = 1, nbpac 
    346             zh_newice(ji)     = hiccrit(1) 
    347          END DO 
    348          IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:) 
     328            zh_newice(ji) = hiccrit(1) 
     329         END DO 
     330         IF( fraz_swi == 1.0 )  zh_newice(:) = hicol_b(:) 
    349331 
    350332         !---------------------- 
    … …  
    352334         !---------------------- 
    353335 
    354          IF ( num_sal .EQ. 1 ) THEN 
    355             zs_newice(:)      =   bulk_sal 
    356          ENDIF ! num_sal 
    357  
    358          IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN 
    359  
    360             DO ji = 1, nbpac 
    361                zs_newice(ji)  =   MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max ) 
    362                zji            =   MOD( npac(ji) - 1, jpi ) + 1 
    363                zjj            =   ( npac(ji) - 1 ) / jpi + 1 
    364                zs_newice(ji)  =   MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) ) 
    365             END DO ! jl 
    366  
    367          ENDIF ! num_sal 
    368  
    369          IF ( num_sal .EQ. 3 ) THEN 
    370             zs_newice(:)      =   2.3 
    371          ENDIF ! num_sal 
     336         SELECT CASE ( num_sal ) 
     337         CASE ( 1 )                    ! Sice = constant  
     338            zs_newice(:) = bulk_sal 
     339         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
     340            DO ji = 1, nbpac 
     341               zji =   MOD( npac(ji) - 1 , jpi ) + 1 
     342               zjj =      ( npac(ji) - 1 ) / jpi + 1 
     343               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj)  ) 
     344            END DO 
     345         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     346            zs_newice(:) =   2.3 
     347         END SELECT 
     348 
    372349 
    373350         !------------------------- 
    … …  
    376353         ! We assume that new ice is formed at the seawater freezing point 
    377354         DO ji = 1, nbpac 
    378             ztmelts           = - tmut * zs_newice(ji) + rtt ! Melting point (K) 
    379             ze_newice(ji)     =   rhoic * ( cpic * ( ztmelts - t_bo_b(ji) )    & 
    380                + lfus * ( 1.0 - ( ztmelts - rtt )   & 
    381                / ( t_bo_b(ji) - rtt ) )           & 
    382                - rcp * ( ztmelts-rtt ) ) 
    383             ze_newice(ji)     =   MAX( ze_newice(ji) , 0.0 ) +                 & 
    384                MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) )   &  
    385                * rhoic * lfus 
     355            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) 
     356            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
     357               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
     358               &                       - rcp  *         ( ztmelts - rtt )  ) 
     359            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    & 
     360               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus 
    386361         END DO ! ji 
    387362         !---------------- 
    … …  
    389364         !---------------- 
    390365         DO ji = 1, nbpac 
    391             zo_newice(ji)     = 0.0 
     366            zo_newice(ji) = 0._wp 
    392367         END DO ! ji 
    393368 
    … …  
    396371         !-------------------------- 
    397372         DO ji = 1, nbpac 
    398             zqbgow(ji)        = qldif_1d(ji) - qcmif_1d(ji) !<0 
     373            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)    !<0 
    399374         END DO ! ji 
    400375 
    … …  
    403378         !------------------- 
    404379         DO ji = 1, nbpac 
    405             zv_newice(ji)     = - zqbgow(ji) / ze_newice(ji) 
     380            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 
    406381 
    407382            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    408             zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) )     &  
    409                + 1.0 ) / 2.0 * maxfrazb 
    410             zdh_frazb(ji) = zfrazb*zv_newice(ji) 
     383            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
     384            zdh_frazb(ji) =         zfrazb   * zv_newice(ji) 
    411385            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    412386         END DO 
    … …  
    415389         ! Salt flux due to new ice growth 
    416390         !--------------------------------- 
    417          IF ( ( num_sal .EQ. 4 ) ) THEN  
    418             DO ji = 1, nbpac 
    419                zji            = MOD( npac(ji) - 1, jpi ) + 1 
    420                zjj            = ( npac(ji) - 1 ) / jpi + 1 
    421                fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    422                   ( sss_m(zji,zjj) - bulk_sal      ) * rhoic *       & 
    423                   zv_newice(ji) / rdt_ice 
    424             END DO 
    425          ELSE 
    426             DO ji = 1, nbpac 
    427                zji            = MOD( npac(ji) - 1, jpi ) + 1 
    428                zjj            = ( npac(ji) - 1 ) / jpi + 1 
    429                fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    430                   ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic *       & 
    431                   zv_newice(ji) / rdt_ice 
    432             END DO ! ji 
    433          ENDIF 
     391         ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 
     392         DO ji = 1, nbpac 
     393            fseqv_1d(ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
     394         END DO ! ji 
    434395 
    435396         !------------------------------------ 
    … …  
    448409 
    449410         ! keep new ice volume in memory 
    450          CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , & 
    451             jpi, jpj ) 
     411         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj ) 
    452412 
    453413         !----------------- 
    … …  
    459419            zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    460420            zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    461             diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 
     421            diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 
    462422         END DO !ji 
    463423 
    … …  
    628588         ! Update salinity 
    629589         !----------------- 
    630          IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     590         IF(  num_sal == 2  ) THEN      ! Sice = F(z,t) 
    631591            DO jl = 1, jpl 
    632592               DO ji = 1, nbpac 
    633                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes 
     593                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )   ! 0 if no ice and 1 if yes 
    634594                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    635595                  zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb 
    … …  
    645605            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 
    646606            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 
    647             IF (  num_sal == 2  .OR.  num_sal == 4  )   & 
     607            IF (  num_sal == 2  )   & 
    648608               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 
    649609            DO jk = 1, nlay_i 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r3294 r3517  
    4444      !! ** Purpose :   computes new salinities in the ice 
    4545      !! 
    46       !! ** Method  :  4 possibilities 
    47       !!               -> num_sal = 1 -> constant salinity for z,t 
    48       !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005] 
    49       !!               -> num_sal = 3 -> S = S(z)   [multiyear ice] 
    50       !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74] 
     46      !! ** Method  :  3 possibilities 
     47      !!               -> num_sal = 1 -> Sice = cst    [ice salinity constant in both time & space]  
     48      !!               -> num_sal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005] 
     49      !!               -> num_sal = 3 -> Sice = S(z)   [multiyear ice] 
    5150      !!--------------------------------------------------------------------- 
    52       INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index 
     51      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index 
    5352      ! 
    5453      INTEGER  ::   ji, jk     ! dummy loop indices  
    55       INTEGER  ::   zji, zjj   ! local integers 
    5654      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars 
    5755      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars 
    … …  
    6462      ! 1) Constant salinity, constant in time                                       | 
    6563      !------------------------------------------------------------------------------| 
    66 !!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !! 
    67       IF( num_sal == 1 ) THEN 
    68          ! 
    69          DO jk = 1, nlay_i 
    70             DO ji = kideb, kiut 
    71                s_i_b(ji,jk) =  bulk_sal 
    72             END DO ! ji 
    73          END DO ! jk 
    74          ! 
    75          DO ji = kideb, kiut 
    76             sm_i_b(ji)      =  bulk_sal  
    77          END DO ! ji 
    78          ! 
     64!!gm comment: if num_sal = 1 s_i_new, s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !! 
     65!!gm           ===>>>   simplification of almost all test on num_sal value 
     66      IF(  num_sal == 1  ) THEN 
     67            s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal 
     68            sm_i_b (kideb:kiut)          =  bulk_sal  
     69            s_i_new(kideb:kiut)          =  bulk_sal 
    7970      ENDIF 
    8071 
    … …  
    8374      !------------------------------------------------------------------------------| 
    8475 
    85       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     76      IF(  num_sal == 2  ) THEN 
    8677 
    8778         !--------------------------------- 
    … …  
    118109            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  
    119110            !                                   ! drainage by flushing   
    120             dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
     111            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
    121112 
    122113            !----------------- 
    … …  
    133124         END DO ! ji 
    134125 
    135          ! Salinity profile 
    136          CALL lim_var_salprof1d( kideb, kiut ) 
     126         CALL lim_var_salprof1d( kideb, kiut )         ! Salinity profile 
     127 
    137128 
    138129         !---------------------------- 
    … …  
    143134!!gm useless 
    144135            ! iflush  : 1 if summer  
    145             iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) )  
     136            iflush  =  MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt        ) )  
    146137            ! igravdr : 1 if t_su lt t_bo 
    147             igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
     138            igravdr =  MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
    148139            ! iaccrbo : 1 if bottom accretion 
    149             iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) ) 
     140            iaccrbo =  MAX( 0._wp , SIGN( 1._wp , dh_i_bott(ji)           ) ) 
    150141!!gm end useless 
    151142            ! 
    … …  
    157148         !---------------------------- 
    158149         DO ji = kideb, kiut 
    159             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
     150            i_ice_switch = 1._wp - MAX(  0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
    160151            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         & 
    161                &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice 
    162             IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 
     152               &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) * r1_rdtice 
    163153         END DO ! ji 
    164154 
    … …  
    179169         END DO 
    180170         ! 
    181       ENDIF ! num_sal .EQ. 2 
     171      ENDIF  
    182172 
    183173      !------------------------------------------------------------------------------| 
    … …  
    185175      !------------------------------------------------------------------------------| 
    186176 
    187       IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut ) 
    188  
    189       !------------------------------------------------------------------------------| 
    190       !  Module 4 : Constant salinity varying in time                                | 
    191       !------------------------------------------------------------------------------| 
    192  
    193       IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974 
    194          ! 
    195          DO ji = kideb, kiut 
    196             zsold = sm_i_b(ji) 
    197             IF( ht_i_b(ji) < 0.4 ) THEN 
    198                sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji)  
    199             ELSE 
    200                sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji) 
    201                sm_i_b(ji) = MIN( sm_i_b(ji) , zsold )   
    202             ENDIF 
    203             IF( ht_i_b(ji) > 3.06918239 ) THEN  
    204                sm_i_b(ji) = 3._wp 
    205             ENDIF 
    206             DO jk = 1, nlay_i 
    207                s_i_b(ji,jk)   = sm_i_b(ji) 
    208             END DO 
    209          END DO 
    210          ! 
    211       ENDIF ! num_sal 
     177      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut ) 
     178 
    212179 
    213180      !------------------------------------------------------------------------------| 
    214181      ! 5) Computation of salt flux due to Bottom growth 
    215182      !------------------------------------------------------------------------------| 
    216  
    217       IF ( num_sal == 4 ) THEN 
    218          DO ji = kideb, kiut 
    219             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    220             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    221             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    )               & 
    222                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice 
    223          END DO 
    224       ELSE 
    225          DO ji = kideb, kiut 
    226             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    227             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    228             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) )               & 
    229                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice 
    230          END DO 
    231       ENDIF 
     183      ! note: s_i_new = bulk_sal in constant salinity case 
     184      DO ji = kideb, kiut 
     185         fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0._wp ) * r1_rdtice 
     186      END DO 
    232187      ! 
    233188      CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold ) 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r3419 r3517  
    392392 
    393393                  ! Ice salinity and age 
    394                   zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , & 
    395                      zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 
    396                   IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    397                      smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0 
    398  
    399                   zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &  
    400                      MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl) 
    401                   oa_i (ji,jj,jl)  = zindic*zage  
     394                  zsal = MAX( MIN(  (rhoic-rhosn)/rhoic*sss_m(ji,jj) ,   & 
     395                     &              zusvoic * zs0sm(ji,jj,jl)         ) , s_i_min ) * v_i(ji,jj,jl) 
     396                  IF(  num_sal == 2  )   smv_i(ji,jj,jl) = zindic * zsal + (1.0-zindic) * 0._wp 
     397 
     398                  zage = MAX(  MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl) 
     399                  oa_i (ji,jj,jl)  = zindic * zage  
    402400 
    403401                  ! Snow heat content 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90

    r3419 r3517  
    190190 
    191191                  ! is there any ice left ? 
    192                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
     192                  zindic = MAX( rzero , SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    193193                  !=1 if hi > 1e-3 and 0 if not 
    194                   zdvres        = MAX(0.0,-v_i(ji,jj,jl)) !residual volume if too much ice was molten 
     194                  zdvres = MAX( 0.0 , -v_i(ji,jj,jl) ) !residual volume if too much ice was molten 
    195195                  !this quantity is positive 
    196                   v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)    !ice volume cannot be negative 
     196                  v_i(ji,jj,jl) = zindic * v_i(ji,jj,jl)    !ice volume cannot be negative 
    197197                  !correct thermodynamic ablation 
    198                   d_v_i_thd(ji,jj,jl)  = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl))  
     198                  d_v_i_thd(ji,jj,jl) = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl))  
    199199                  ! THIS IS NEW 
    200                   d_a_i_thd(ji,jj,jl)  = zindic  *  d_a_i_thd(ji,jj,jl) + &  
    201                      (1.0-zindic) * (-old_a_i(ji,jj,jl))  
     200                  d_a_i_thd(ji,jj,jl) = zindic  *  d_a_i_thd(ji,jj,jl) + (1.0-zindic) * (-old_a_i(ji,jj,jl))  
    202201 
    203202                  !residual salt flux if ice is over-molten 
    204                   fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
    205                      ( rhoic * zdvres / rdt_ice ) 
    206                   !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres / rdt_ice 
     203                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres * r1_rdtice ) 
     204                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres * r1_rdtice 
    207205 
    208206                  ! is there any snow left ? 
    209                   zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )  
    210                   zvsold        = v_s(ji,jj,jl) 
    211                   zdvres        = MAX(0.0,-v_s(ji,jj,jl)) !residual volume if too much ice was molten 
     207                  zindsn = MAX(  rzero, SIGN( rone , v_s(ji,jj,jl) - epsi10 ) )  
     208                  zvsold = v_s(ji,jj,jl) 
     209                  zdvres = MAX( 0.0 , -v_s(ji,jj,jl) )  !residual volume if too much ice was molten 
    212210                  !this quantity is positive 
    213211                  v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl)    !snow volume cannot be negative 
    … …  
    215213                  d_v_s_thd(ji,jj,jl)  = zindsn  *  d_v_s_thd(ji,jj,jl) + &  
    216214                     (1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl))  
    217                   !unsure correction on salt flux.... maybe future will tell it was not that right 
    218  
    219                   !residual salt flux if snow is over-molten 
    220                   fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sss_m(ji,jj) * ( rhosn * zdvres / rdt_ice ) 
    221                   !this flux will be positive if snow was over-molten 
    222                   !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres / rdt_ice 
     215 
     216                  !no salt flux when snow is over-molten 
     217                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres * r1_rdtice 
    223218               ENDIF 
    224219            END DO !ji 
    … …  
    306301      !-------------- 
    307302 
    308       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case 
     303      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)] 
    309304         ! 
    310305         IF( ln_nicep ) THEN   
    … …  
    317312            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl) 
    318313         ENDIF 
    319  
     314         ! 
    320315         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) + d_smv_i_trp(:,:,:) 
    321316         ! 
    … …  
    616611            DO ji = 1, jpi 
    617612               IF ( internal_melt(ji,jj,jl) == 1 ) THEN 
    618                   v_s(ji,jj,jl)   = 0.0 
    619                   e_s(ji,jj,1,jl) = 0.0 
     613                  v_s(ji,jj,jl)   = 0._wp 
     614                  e_s(ji,jj,1,jl) = 0._wp 
    620615                  !   ! release heat 
    621                   fheat_res(ji,jj) = fheat_res(ji,jj)  & 
    622                      + ze_s * v_s(ji,jj,jl) / rdt_ice 
     616                  fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s  * v_s(ji,jj,jl) * r1_rdtice 
    623617                  ! release mass 
    624                   rdm_snw(ji,jj) =  rdm_snw(ji,jj) + rhosn * v_s(ji,jj,jl) 
     618                  rdm_snw  (ji,jj) = rdm_snw  (ji,jj) + rhosn * v_s(ji,jj,jl) 
    625619               ENDIF 
    626620            END DO 
    … …  
    648642               !                ENDIF 
    649643               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 
    650                   oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl) 
     644                  oa_i(ji,jj,jl) = rdt_ice * numit / 86400.0 * a_i(ji,jj,jl) 
    651645               ENDIF 
    652646               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 
    … …  
    657651               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)  
    658652               ! snow thickness cannot be smaller than 1e-6 
    659                v_s(ji,jj,jl)  = zindsn*v_s(ji,jj,jl)*zindb 
     653               v_s(ji,jj,jl)  = zindsn * v_s(ji,jj,jl) * zindb 
    660654               v_s(ji,jj,jl)  = v_s(ji,jj,jl) *  MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) ) 
    661655 
    … …  
    737731      !--------------------- 
    738732 
    739       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case 
    740  
     733      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)] 
     734         ! 
    741735         DO jl = 1, jpl 
    742736            DO jk = 1, nlay_i 
    743737               DO jj = 1, jpj  
    744                   DO ji = 1, jpi 
    745                      ! salinity stays in bounds 
    746                      smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), & 
    747                         0.1 * v_i(ji,jj,jl) ) 
     738                  DO ji = 1, jpi           ! salinity stays in bounds 
     739                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), 0.1 * v_i(ji,jj,jl) ) 
    748740                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) 
    749                      smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + & 
    750                         0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 
     741                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 
    751742                  END DO ! ji 
    752743               END DO ! jj 
    753744            END DO !jk 
    754745         END DO !jl 
    755  
     746         ! 
    756747      ENDIF 
    757748 
    … …  
    819810            ! 1) Count the number of existing categories 
    820811            DO jl  = 1, jpl 
     812!!cr : comment the second line of zindb definition, and use epsi04 in the 1st one 
    821813               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) )  
    822814               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) )  

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r3294 r3517  
    182182      END DO 
    183183 
    184       IF(  num_sal == 2  .OR.  num_sal == 4  )THEN 
     184      IF(  num_sal == 2  )THEN 
    185185         DO jl = 1, jpl 
    186186            DO jj = 1, jpj 
    … …  
    309309      ! Vertically constant, constant in time 
    310310      !--------------------------------------- 
    311       IF( num_sal == 1 )   s_i(:,:,:,:) = bulk_sal 
     311      IF(  num_sal == 1 )   s_i(:,:,:,:) = bulk_sal 
    312312 
    313313      !----------------------------------- 
    314314      ! Salinity profile, varying in time 
    315315      !----------------------------------- 
    316  
    317       IF(   num_sal == 2  .OR.   num_sal == 4   ) THEN 
     316      IF(  num_sal == 2  ) THEN 
    318317         ! 
    319318         DO jk = 1, nlay_i 
    … …  
    331330         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf 
    332331         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 
    333  
     332         ! 
    334333         zalpha(:,:,:) = 0._wp 
    335334         DO jl = 1, jpl 
    … …  
    347346            END DO 
    348347         END DO 
    349  
     348         ! 
    350349         dummy_fac = 1._wp / nlay_i                   ! Computation of the profile 
    351350         DO jl = 1, jpl 
    … …  
    361360            END DO ! jk 
    362361         END DO ! jl 
    363  
     362         ! 
    364363      ENDIF ! num_sal 
    365364 
    … …  
    368367      !------------------------------------------------------- 
    369368 
    370       IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
     369      IF(  num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
    371370         ! 
    372371         sm_i(:,:,:) = 2.30_wp 
    … …  
    380379            END DO 
    381380         END DO 
    382  
     381         ! 
    383382      ENDIF ! num_sal 
    384383      ! 
    … …  
    447446      !------------------------------------------------------ 
    448447 
    449       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     448      IF(  num_sal == 2  ) THEN 
    450449         ! 
    451450         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90

    r3402 r3517  
    111111         zcmo(ji,jj,13) = qns(ji,jj) 
    112112         ! See thersf for the coefficient 
    113          zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     113         zcmo(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in Kg/m2/day = mm/day 
    114114         zcmo(ji,jj,15) = utau_ice(ji,jj) 
    115115         zcmo(ji,jj,16) = vtau_ice(ji,jj) 
    … …  
    154154               rcmoy(ji,jj,13) = qns(ji,jj) 
    155155               ! See thersf for the coefficient 
    156                rcmoy(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     156               rcmoy(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in mm/day 
    157157               rcmoy(ji,jj,15) = utau_ice(ji,jj) 
    158158               rcmoy(ji,jj,16) = vtau_ice(ji,jj) 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r2715 r3517  
    6666   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld 
    6767   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif 
    68    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmicif_1d    !: <==> the 2D  rdmicif 
    69    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmsnif_1d    !: <==> the 2D  rdmsnif 
     68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_ice_1d    !: <==> the 2D  rdm_ice 
     69   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_snw_1d    !: <==> the 2D  rdm_snw 
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq 
    7171   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi 
    … …  
    160160      ! 
    161161      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     & 
    162          &      fbif_1d    (jpij) , rdmicif_1d (jpij) , rdmsnif_1d (jpij) ,     & 
     162         &      fbif_1d    (jpij) , rdm_ice_1d (jpij) , rdm_snw_1d (jpij) ,     & 
    163163         &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     & 
    164164         &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     & 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r3402 r3517  
    175175         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate arrays' ) 
    176176         ! 
    177          sfx(:,:) = 0.0_wp                      ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
     177         sfx(:,:) = 0._wp                       ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    178178         ! 
    179179      ENDIF 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r3489 r3517  
    179179         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 
    180180         ! 
    181          sfx(:,:) = 0.0_wp                         ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
     181         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    182182         ! 
    183183      ENDIF 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r3488 r3517  
    779779         ! Stress module can be negative when received (interpolation problem) 
    780780         IF( llnewtau ) THEN  
    781             frcv(jpr_taum)%z3(:,:,1) = MAX( 0.0e0, frcv(jpr_taum)%z3(:,:,1) ) 
     781            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) ) 
    782782         ENDIF 
    783783      ENDIF 
    … …  
    823823         !                                                   ! ========================= ! 
    824824         ! 
    825          !                                                       ! total freshwater fluxes over the ocean (emp, sfx ) 
     825         !                                                       ! total freshwater fluxes over the ocean (emp) 
    826826         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation 
    827827         CASE( 'conservative' ) 
    … …  
    12611261      !                                                      ! ========================= ! 
    12621262      CASE( 'oce only' ) 
    1263          qsr_tot(:,:  ) = MAX(0.0,frcv(jpr_qsroce)%z3(:,:,1)) 
     1263         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) ) 
    12641264      CASE( 'conservative' ) 
    12651265         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
    … …  
    13641364            ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
    13651365         CASE( 'no' ) 
    1366             ztmp3(:,:,:) = 0.0 
     1366            ztmp3(:,:,:) = 0._wp 
    13671367            DO jl=1,jpl 
    13681368               ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl) 
    … …  
    14161416            ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl) 
    14171417         CASE( 'no' ) 
    1418             ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0 
     1418            ztmp3(:,:,:) = 0._wp   ;  ztmp4(:,:,:) = 0._wp 
    14191419            DO jl=1,jpl 
    14201420               ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl) 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r3402 r3517  
    1010   !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine 
    1111   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step 
    12    !!            4.0  ! 2011-01  (A Porter)  dynamical allocation 
     12   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation 
    1313   !!---------------------------------------------------------------------- 
    1414#if defined key_lim3 
    … …  
    170170 
    171171         !                                           ! intialisation to zero    !!gm is it truly necessary ??? 
    172          d_a_i_thd  (:,:,:)   = 0.e0   ;   d_a_i_trp  (:,:,:)   = 0.e0 
    173          d_v_i_thd  (:,:,:)   = 0.e0   ;   d_v_i_trp  (:,:,:)   = 0.e0 
    174          d_e_i_thd  (:,:,:,:) = 0.e0   ;   d_e_i_trp  (:,:,:,:) = 0.e0 
    175          d_v_s_thd  (:,:,:)   = 0.e0   ;   d_v_s_trp  (:,:,:)   = 0.e0 
    176          d_e_s_thd  (:,:,:,:) = 0.e0   ;   d_e_s_trp  (:,:,:,:) = 0.e0 
    177          d_smv_i_thd(:,:,:)   = 0.e0   ;   d_smv_i_trp(:,:,:)   = 0.e0 
    178          d_oa_i_thd (:,:,:)   = 0.e0   ;   d_oa_i_trp (:,:,:)   = 0.e0 
    179          ! 
    180          fseqv    (:,:) = 0.e0 
    181          fsbri    (:,:) = 0.e0     ;   fsalt_res(:,:) = 0.e0 
    182          fsalt_rpo(:,:) = 0.e0 
    183          fhmec    (:,:) = 0.e0     ;   fhbri    (:,:) = 0.e0 
    184          fmmec    (:,:) = 0.e0     ;   fheat_res(:,:) = 0.e0 
    185          fheat_rpo(:,:) = 0.e0     ;   focea2D  (:,:) = 0.e0 
    186          fsup2D   (:,:) = 0.e0 
     172         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
     173         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp 
     174         d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp 
     175         d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp 
     176         d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp 
     177         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
     178         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
     179         ! 
     180         fseqv    (:,:) = 0._wp 
     181         fsbri    (:,:) = 0._wp     ;   fsalt_res(:,:) = 0._wp 
     182         fsalt_rpo(:,:) = 0._wp 
     183         fhmec    (:,:) = 0._wp     ;   fhbri    (:,:) = 0._wp 
     184         fmmec    (:,:) = 0._wp     ;   fheat_res(:,:) = 0._wp 
     185         fheat_rpo(:,:) = 0._wp     ;   focea2D  (:,:) = 0._wp 
     186         fsup2D   (:,:) = 0._wp 
    187187         !  
    188          diag_sni_gr(:,:) = 0.e0   ;   diag_lat_gr(:,:) = 0.e0 
    189          diag_bot_gr(:,:) = 0.e0   ;   diag_dyn_gr(:,:) = 0.e0 
    190          diag_bot_me(:,:) = 0.e0   ;   diag_sur_me(:,:) = 0.e0 
     188         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp 
     189         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp 
     190         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp 
    191191         ! dynamical invariants 
    192          delta_i(:,:) = 0.e0       ;   divu_i(:,:) = 0.e0       ;   shear_i(:,:) = 0.e0 
     192         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp 
    193193 
    194194                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    … …  
    196196         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    197197         ! 
    198          IF( .NOT. lk_c1d ) THEN 
    199                                                      ! Ice dynamics & transport (not in 1D case) 
     198         IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case) 
    200199                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    201200                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    … …  
    210209                          CALL lim_var_bv                 ! bulk brine volume (diag) 
    211210                          CALL lim_thd( kt )              ! Ice thermodynamics  
    212                           zcoef = rdt_ice / 86400.e0      !  Ice natural aging 
     211                          zcoef = rdt_ice /rday           !  Ice natural aging 
    213212                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    214213                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    … …  
    268267 
    269268      inb_altests = 10 
    270       inb_alp(:)  = 0 
     269      inb_alp(:)  =  0 
    271270 
    272271      ! Alert if incompatible volume and concentration 
    … …  
    277276         DO jj = 1, jpj 
    278277            DO ji = 1, jpi 
    279                IF(  v_i(ji,jj,jl) /= 0.e0   .AND.   a_i(ji,jj,jl) == 0.e0   ) THEN 
     278               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    280279                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    281280                  WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    … …  
    297296      DO jj = 1, jpj 
    298297         DO ji = 1, jpi 
    299             IF(   ht_i(ji,jj,jl) .GT. 50.0   ) THEN 
     298            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    300299               CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    301300               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    309308      DO jj = 1, jpj 
    310309         DO ji = 1, jpi 
    311             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) .GT. 0.5  .AND.  & 
    312                &  at_i(ji,jj) .GT. 0.e0   ) THEN 
     310            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  & 
     311               &  at_i(ji,jj) > 0._wp   ) THEN 
    313312               CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    314313               WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    … …  
    332331      DO jj = 1, jpj 
    333332         DO ji = 1, jpi 
    334             IF(   tms(ji,jj) .LE. 0.0   .AND.   at_i(ji,jj) .GT. 0.e0   ) THEN  
     333            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    335334               CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    336335               WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    … …  
    356355            DO ji = 1, jpi 
    357356!!gm  test twice sm_i ...  ????  bug? 
    358                IF( ( ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) .OR. & 
    359                      ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) ) .AND. & 
    360                              ( a_i(ji,jj,jl) .GT. 0.e0 ) ) THEN 
     357               IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR. & 
     358                     ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. & 
     359                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    361360!                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    362361!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    … …  
    377376         DO jj = 1, jpj 
    378377            DO ji = 1, jpi 
    379                IF ( ( ( ABS( o_i(ji,jj,jl) ) .GT. rdt_ice ) .OR. & 
    380                       ( ABS( o_i(ji,jj,jl) ) .LT. 0.00) ) .AND. & 
    381                              ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN 
     378               IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. & 
     379                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
     380                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    382381                  CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    383382                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    412411      DO jj = 1, jpj 
    413412         DO ji = 1, jpi 
    414             IF(   ABS( qns(ji,jj) ) .GT. 1500.0   .AND.  ( at_i(ji,jj) .GT. 0.0 ) )  THEN 
     413            IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp )  THEN 
    415414               ! 
    416415               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    … …  
    450449               DO ji = 1, jpi 
    451450                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt 
    452                   IF( t_i(ji,jj,jk,jl) .GE. ztmelts  .AND.  v_i(ji,jj,jl) .GT. 1.e-6   & 
    453                      &                               .AND.  a_i(ji,jj,jl) .GT. 0.e0    ) THEN 
     451                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   & 
     452                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    454453                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    455454                     WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 

  • branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r3396 r3517  
    8282   END FUNCTION sbc_rnf_alloc 
    8383 
     84 
    8485   SUBROUTINE sbc_rnf( kt ) 
    8586      !!---------------------------------------------------------------------- 
    … …  
    9596      !!---------------------------------------------------------------------- 
    9697      INTEGER, INTENT(in) ::   kt          ! ocean time step 
    97       !! 
     98      ! 
    9899      INTEGER  ::   ji, jj   ! dummy loop indices 
    99100      !!---------------------------------------------------------------------- 
    … …  
    126127         ! 
    127128         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    128             rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )   
     129            ! 
     130            rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt 
    129131            ! 
    130132            !                                                     ! set temperature & salinity content of runoffs 
    … …  
    248250      INTEGER           ::   ji, jj, jk    ! dummy loop indices 
    249251      INTEGER           ::   ierror, inum  ! temporary integer 
    250       !!  
     252      ! 
    251253      NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   & 
    252254         &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &   
    253255         &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact   
    254256      !!---------------------------------------------------------------------- 
    255  
     257      ! 
    256258      !                                   ! ============ 
    257259      !                                   !   Namelist 
    … …  
    269271      REWIND ( numnam )                         ! Read Namelist namsbc_rnf 
    270272      READ   ( numnam, namsbc_rnf ) 
    271  
     273      ! 
    272274      !                                         ! Control print 
    273275      IF(lwp) THEN 
    … …  
    282284         WRITE(numout,*) '      multiplicative factor for runoff           rn_rfact     = ', rn_rfact     
    283285      ENDIF 
    284  
     286      ! 
    285287      !                                   ! ================== 
    286288      !                                   !   Type of runoff 
    … …  
    391393            nkrnf = 2 
    392394            DO WHILE( nkrnf /= jpkm1 .AND. gdepw_0(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1   ;   END DO 
    393             IF( ln_sco )   & 
    394                CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' ) 
     395            IF( ln_sco )   CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' ) 
    395396         ENDIF 
    396397         IF(lwp) WRITE(numout,*) 
    … …  
    410411         nkrnf = 0 
    411412      ENDIF 
    412  
     413      ! 
    413414   END SUBROUTINE sbc_rnf_init 
    414415 
    … …  
    434435      !!                rnfmsk_z vertical structure 
    435436      !!---------------------------------------------------------------------- 
    436       ! 
    437437      INTEGER            ::   inum        ! temporary integers 
    438438      CHARACTER(len=140) ::   cl_rnfile   ! runoff file name 
    … …  
    442442      IF(lwp) WRITE(numout,*) 'rnf_mouth : river mouth mask' 
    443443      IF(lwp) WRITE(numout,*) '~~~~~~~~~ ' 
    444  
     444      ! 
    445445      cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname ) 
    446446      IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4)' ) TRIM( cl_rnfile ), nyear    ! add year 
    447447         IF( sn_cnf%cltype == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month 
    448448      ENDIF 
    449    
     449      ! 
    450450      ! horizontal mask (read in NetCDF file) 
    451451      CALL iom_open ( cl_rnfile, inum )                           ! open file 
    452452      CALL iom_get  ( inum, jpdom_data, sn_cnf%clvar, rnfmsk )    ! read the river mouth array 
    453453      CALL iom_close( inum )                                      ! close file 
    454        
     454      ! 
    455455      IF( nclosea == 1 )    CALL clo_rnf( rnfmsk )                ! closed sea inflow set as ruver mouth 
    456  
     456      ! 
    457457      rnfmsk_z(:)   = 0._wp                                        ! vertical structure  
    458458      rnfmsk_z(1)   = 1.0 
Note: See TracChangeset for help on using the changeset viewer.