Changeset 4220


Ignore:
Timestamp:
2013-11-15T16:36:52+01:00 (8 years ago)
Author:
clem
Message:

corrections for restartability

Location:
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4155 r4220  
    399399   REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean 
    400400   REAL(wp)              , PUBLIC ::   amax          = 0.99               !: maximum ice concentration 
    401    !                                                                      ! 
     401   ! 
    402402   !!-------------------------------------------------------------------------- 
    403403   !! * Ice diagnostics 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r4155 r4220  
    8787      CALL lim_itd_ini                 ! ice thickness distribution initialization 
    8888      ! 
    89       CALL lim_sbc_init                ! ice surface boundary condition    
    90  
    91  
    9289      !                                ! Initial sea-ice state 
    9390      IF( .NOT.ln_rstart ) THEN              ! start from rest 
     
    104101      ENDIF 
    105102      ! 
     103      CALL lim_sbc_init                ! ice surface boundary condition    
     104      ! 
    106105      fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction 
     106      tn_ice(:,:,:) = t_su(:,:,:) 
    107107      ! 
    108108      nstart = numit  + nn_fsbc       
     
    239239      END DO 
    240240      ! 
    241       tn_ice(:,:,:) = t_su(:,:,:) 
    242241      ! 
    243242   END SUBROUTINE lim_itd_ini 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r4155 r4220  
    123123 
    124124      REAL(wp) ::   dtevp              ! time step for subcycling 
    125       REAL(wp) ::   dtotel, ecc2      ! square of yield ellipse eccenticity 
     125      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    126126      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars 
    127127      REAL(wp) ::   zu_ice2, zv_ice1   ! 
     
    199199            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 
    200200#if defined key_lim3 
    201             zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp 
     201            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) 
    202202#endif 
    203203#if defined key_lim2 
     
    325325      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    326326      ecc2 = ecc * ecc 
     327      ecci = 1. / ecc2 
    327328 
    328329      !-Initialise stress tensor  
     
    425426               !-Calculate stress tensor components zs1 and zs2  
    426427               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    427                zs1(ji,jj) = ( zs1(ji,jj) & 
    428                   &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    & 
    429                   &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) & 
    430                   * zpresh(ji,jj) ) )                          &        
    431                   &        / ( 1.0 + alphaevp * dtotel ) 
    432  
    433                zs2(ji,jj) = ( zs2(ji,jj)   & 
    434                   &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  & 
    435                   zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) & 
    436                   &        / ( 1.0 + alphaevp*ecc2*dtotel ) 
     428               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
     429                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
     430                  &          / ( 1._wp + alphaevp * dtotel ) 
     431 
     432               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
     433                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
     434                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
     435 
     436               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
     437               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
     438               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
     439               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
     440               !   &         / ( 1._wp + dtotel ) 
    437441 
    438442            END DO 
     
    468472 
    469473               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    470                zs12(ji,jj) = ( zs12(ji,jj)      & 
    471                   &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / & 
    472                   &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) & 
    473                   &         / ( 1.0 + alphaevp*ecc2*dtotel )  
     474               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
     475                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     476                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
     477 
     478               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
     479               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
     480               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     481               !   &          / ( 1.0 + dtotel )  
    474482 
    475483            END DO ! ji 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r4099 r4220  
    1717   !!---------------------------------------------------------------------- 
    1818   USE ice            ! sea-ice variables 
     19   USE oce     , ONLY :  snwice_mass, snwice_mass_b 
    1920   USE par_ice        ! sea-ice parameters 
    2021   USE dom_oce        ! ocean domain 
     
    159160      END DO 
    160161 
    161       CALL iom_rstput( iter, nitrst, numriw, 'u_ice'     , u_ice      ) 
    162       CALL iom_rstput( iter, nitrst, numriw, 'v_ice'     , v_ice      ) 
    163       CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'     , fsbbq      ) 
    164       CALL iom_rstput( iter, nitrst, numriw, 'iatte'     , iatte      ) ! clem modif 
    165       CALL iom_rstput( iter, nitrst, numriw, 'oatte'     , oatte      ) ! clem modif 
    166       CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 
    167       CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  ) 
    168       CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 
     162      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'        , u_ice      ) 
     163      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'        , v_ice      ) 
     164      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'        , fsbbq      ) 
     165      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'    , stress1_i  ) 
     166      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  ) 
     167      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i ) 
     168      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass )   !clem modif 
     169      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) !clem modif 
    169170 
    170171      DO jl = 1, jpl  
     
    371372         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    372373         t_su(:,:,jl) = z2d(:,:) 
    373          tn_ice (:,:,:) = t_su (:,:,:) 
    374       END DO 
    375  
    376       DO jl = 1, jpl  
    377          CALL lbc_lnk( smv_i(:,:,jl) , 'T' ,  1. ) 
    378          CALL lbc_lnk( v_i  (:,:,jl) , 'T' ,  1. ) 
    379          CALL lbc_lnk( a_i  (:,:,jl) , 'T' ,  1. ) 
    380       END DO 
    381  
    382       ! we first with bulk ice salinity 
    383       DO jl = 1, jpl 
    384          DO jj = 1, jpj 
    385             DO ji = 1, jpi 
    386                zindb          = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - 1.0e-4 ) )  
    387                sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),1.0e-6) * zindb 
    388                ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX(a_i(ji,jj,jl),1.0e-6) * zindb 
    389             END DO 
    390          END DO 
    391       END DO 
    392  
    393       DO jk = 1, nlay_i 
    394          s_i(:,:,jk,:) = sm_i(:,:,:) 
    395       END DO 
    396  
    397       IF( num_sal == 2 ) THEN      ! Salinity profile 
    398          DO jl = 1, jpl 
    399             DO jk = 1, nlay_i 
    400                DO jj = 1, jpj 
    401                   DO ji = 1, jpi 
    402                      zs_inf        = sm_i(ji,jj,jl) 
    403                      z_slope_s     = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01_wp , ht_i(ji,jj,jl) ) 
    404                      !- slope of the salinity profile 
    405                      zs_zero(jk)   = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) / REAL(nlay_i,wp) 
    406                      zsmax = 4.5_wp 
    407                      zsmin = 3.5_wp 
    408                      IF(     sm_i(ji,jj,jl) < zsmin ) THEN 
    409                         zalpha = 1._wp 
    410                      ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN 
    411                         zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin ) 
    412                      ELSE 
    413                         zalpha = 0._wp 
    414                      ENDIF 
    415                      s_i(ji,jj,jk,jl) = zalpha * zs_zero(jk) + ( 1._wp - zalpha ) * zs_inf 
    416                   END DO 
    417                END DO 
    418             END DO 
    419          END DO 
    420       ENDIF 
     374      END DO 
    421375 
    422376      DO jl = 1, jpl  
     
    440394      CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      ) 
    441395      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'     , fsbbq      ) 
    442       CALL iom_get( numrir, jpdom_autoglo, 'iatte'     , iatte      ) ! clem modif 
    443       CALL iom_get( numrir, jpdom_autoglo, 'oatte'     , oatte      ) ! clem modif 
    444396      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  ) 
    445397      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
    446398      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 
     399      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass )   !clem modif 
     400      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) !clem modif 
    447401 
    448402      DO jl = 1, jpl  
     
    568522      END DO 
    569523      ! 
    570       !clem CALL iom_close( numrir ) 
     524      !CALL iom_close( numrir ) !clem: closed in sbcice_lim.F90 
    571525      ! 
    572526      CALL wrk_dealloc( nlay_i, zs_zero ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4099 r4220  
    3838   USE cpl_oasis3, ONLY : lk_cpl 
    3939   USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget 
    40    USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 
     40   USE oce,        ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 
    4141   USE dom_ice,    ONLY : tms 
    4242   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     
    249249            ! mass flux at the ocean/ice interface (sea ice fraction) 
    250250            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean 
    251             zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting   
     251            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus melting   
    252252 
    253253            fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model 
     
    416416      ENDIF 
    417417      ! clem modif 
    418       iatte(:,:) = 1._wp 
    419       oatte(:,:) = 1._wp 
    420       ! 
    421       !                                      ! embedded sea ice 
    422       IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
    423          snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
    424          snwice_mass_b(:,:) = snwice_mass(:,:) 
    425       ELSE 
    426          snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges 
    427          snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges 
    428       ENDIF 
    429       IF( nn_ice_embd == 2  .AND.         &  ! full embedment (case 2) & no restart 
    430          &  .NOT. ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area 
    431          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    432          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
    433          ! 
    434          ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n]  
    435          !       which were previously set in domvvl 
    436          IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only??? 
    437             DO jj = 1, jpjm1 
    438                DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
    439                   zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    440                   zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    441                   zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
    442                   sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    443                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    444                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    445                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    446                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    447                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    448                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    449                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     418      IF( .NOT. ln_rstart ) THEN 
     419         iatte(:,:) = 1._wp 
     420         oatte(:,:) = 1._wp 
     421      ENDIF 
     422      ! 
     423      ! clem: snwice_mass in the restart file now 
     424      IF( .NOT. ln_rstart ) THEN 
     425         !                                      ! embedded sea ice 
     426         IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
     427            snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
     428            snwice_mass_b(:,:) = snwice_mass(:,:) 
     429         ELSE 
     430            snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges 
     431            snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges 
     432         ENDIF 
     433         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area 
     434            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
     435            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
     436            ! 
     437            ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n]  
     438            !       which were previously set in domvvl 
     439            IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only??? 
     440               DO jj = 1, jpjm1 
     441                  DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
     442                     zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     443                     zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     444                     zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     445                     sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     446                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
     447                     sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     448                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
     449                     sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     450                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
     451                     sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
     452                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     453                  END DO 
    450454               END DO 
    451             END DO 
    452             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
    453             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
    454             DO jj = 1, jpjm1 
    455                DO ji = 1, jpim1      ! NO Vector Opt. 
    456                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    457                        &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    458                        &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    459                        &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
     455               CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     456               CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     457               DO jj = 1, jpjm1 
     458                  DO ji = 1, jpim1      ! NO Vector Opt. 
     459                     sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     460                          &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     461                          &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
     462                          &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
     463                  END DO 
    460464               END DO 
    461             END DO 
    462             CALL lbc_lnk( sshf_n, 'F', 1. ) 
    463           ENDIF 
    464       ENDIF 
     465               CALL lbc_lnk( sshf_n, 'F', 1. ) 
     466            ENDIF 
     467         ENDIF 
     468      ENDIF ! .NOT. ln_rstart 
    465469      ! 
    466470!!?      IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4072 r4220  
    2222   USE phycst         ! physical constants 
    2323   USE dom_oce        ! ocean space and time domain variables 
     24   USE oce     , ONLY :  iatte, oatte 
    2425   USE ice            ! LIM: sea-ice variables 
    2526   USE par_ice        ! LIM: sea-ice parameters 
     
    220221            ! 
    221222            ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 
    222             qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 
     223            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 
    223224            ! 
    224225            ! oceanic heat flux (limthd_dh) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4045 r4220  
    656656         DO ji = kideb , kiut 
    657657            ! snow temperatures       
    658             IF (ht_s_b(ji).GT.0) & 
     658            IF (ht_s_b(ji).GT.0._wp) & 
    659659               t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    660660               *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 
     
    692692            DO ji = kideb , kiut 
    693693               ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt  
    694                t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0) 
     694               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp) 
    695695               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) 
    696696            END DO 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r4155 r4220  
    7474      REAL(wp) ::   zindb  , zindsn , zindic, zindh, zinda      ! local scalar 
    7575      REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      - 
    76       REAL(wp) ::   zcfl , zusnit , zrtt          !   -      - 
     76      REAL(wp) ::   zcfl , zusnit                 !   -      - 
    7777      REAL(wp) ::   ze   , zsal   , zage          !   -      - 
    7878      ! 
     
    450450                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 ) 
    451451                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 
    452                   zrtt            = 173.15 * rone  
    453452                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 
    454453                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r4045 r4220  
    180180               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
    181181               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
    182                a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction 
    183182            END DO 
    184183         END DO 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4072 r4220  
    2626   USE restart         ! ocean restart 
    2727   USE wrk_nemo         ! work arrays 
     28   USE sbcrnf         ! river runoffd 
    2829 
    2930   IMPLICIT NONE 
     
    7172      REAL(wp)   ::   z_hc        , z_sc          ! heat and salt content 
    7273      REAL(wp)   ::   z_v1        , z_v2          ! volume 
    73       REAL(wp)   ::   z1_rau0                     ! local scalars 
    7474      REAL(wp)   ::   zdeltat                     !    -     - 
    7575      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
     
    8686      ! 1 - Trends due to forcing ! 
    8787      ! ------------------------- ! 
    88       z1_rau0 = 1.e0 / rau0 
    89       z_frc_trd_v = z1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes 
     88      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes 
    9089      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * zsurf(:,:) )       ! heat fluxes 
    9190      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * zsurf(:,:) )       ! salt fluxes 
     91      ! 
     92      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * zsurf(:,:) ) 
     93      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * zsurf(:,:) ) 
     94 
    9295      ! Add penetrative solar radiation 
    9396      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * zsurf(:,:) ) 
    9497      ! Add geothermal heat flux 
    95       IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qgh_trd0(:,:) * zsurf(:,:) ) 
     98      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * zsurf(:,:) ) 
    9699      ! 
    97100      frc_v = frc_v + z_frc_trd_v * rdt 
     
    124127       
    125128      ! add ssh if not vvl 
    126 #if ! defined key_vvl 
    127      zdiff_v2 = zdiff_v2 + zdiff_v1 
    128      zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem)   & 
    129             &                           - hcssh_loc_ini(:,:) ) ) 
    130      zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal)   & 
    131             &                           - scssh_loc_ini(:,:) ) ) 
    132 #endif  
     129      IF( .NOT. lk_vvl ) THEN 
     130        zdiff_v2 = zdiff_v2 + zdiff_v1 
     131        zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem)   & 
     132               &                           - hcssh_loc_ini(:,:) ) ) 
     133        zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal)   & 
     134               &                           - scssh_loc_ini(:,:) ) ) 
     135      ENDIF 
    133136      ! 
    134137      ! ----------------------- ! 
     
    149152      ENDDO 
    150153      ! add ssh if not vvl 
    151 #if ! defined key_vvl 
    152      z_v2 = z_v2 + z_v1 
    153      z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) ) 
    154      z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) ) 
    155 #endif  
     154      IF( .NOT. lk_vvl ) THEN 
     155        z_v2 = z_v2 + z_v1 
     156        z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) ) 
     157        z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) ) 
     158      ENDIF 
    156159 
    157160      ! ----------------------- ! 
     
    160163      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt ) 
    161164! 
    162       CALL iom_put( 'bgtemper',z_hc / z_v2 )               ! Temperature (C)  
    163       CALL iom_put( 'bgsaline',z_sc / z_v2 )               ! Salinity (psu) 
    164       !CALL iom_put( 'bgheatco',zdiff_hc*fact1*zdeltat )      ! Equivalent heat flux (W/m2) 
    165       !CALL iom_put( 'bgsaltco',zdiff_sc*fact21*zdeltat )     ! Equivalent water flux (mm/s) 
    166       CALL iom_put( 'bgheatco',zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J) 
    167       CALL iom_put( 'bgsaltco',zdiff_sc * 1.e-9 )                        ! Salt content variation (psu*km3)  
    168       CALL iom_put( 'bgvolssh',zdiff_v1 * 1.e-9 )                         ! volume ssh (km3)   
    169       CALL iom_put( 'bgsshtot',zdiff_v1 / glob_sum(zsurf) )              ! ssh (m)   
    170       CALL iom_put( 'bgvoltot',zdiff_v2 * 1.e-9 )                         ! volume total (km3)  
    171       CALL iom_put( 'bgfrcvol',frc_v * 1.e-9 )                         ! vol - surface forcing (volume)  
    172       CALL iom_put( 'bgfrctem',frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content)  
    173       CALL iom_put( 'bgfrcsal',frc_s * 1.e-9 )                         ! sc  - surface forcing (salt content)  
     165      CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C)  
     166      CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu) 
     167      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J) 
     168      CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 )                 ! Salt content variation (psu*km3)  
     169      CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )                    ! volume ssh (km3)   
     170      CALL iom_put( 'bgsshtot' , zdiff_v1 / glob_sum(zsurf) )          ! ssh (m)   
     171      CALL iom_put( 'bgvoltot' , zdiff_v2 * 1.e-9 )                 ! volume total (km3)  
     172      CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 )                     ! vol - surface forcing (volume)  
     173      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content)  
     174      CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content)  
    174175      ! 
    175176      CALL wrk_dealloc( jpi, jpj, zsurf ) 
     
    246247            CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )          
    247248         ENDIF 
    248           
    249          ! ---------------------------------- ! 
    250          ! 4 - initial conservation variables ! 
    251          ! ---------------------------------- ! 
    252          !ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh 
    253          !DO jk = 1, jpk 
    254          !   e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors 
    255          !   hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content 
    256          !   sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content 
    257          !END DO 
    258          !hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
    259          !scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
    260          !frc_v = 0._wp                                           ! volume       trend due to forcing 
    261          !frc_t = 0._wp                                           ! heat content   -    -   -    -    
    262          !frc_s = 0._wp                                           ! salt content   -    -   -    -          
    263249         ! 
    264250         CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r4045 r4220  
    164164 
    165165#endif 
     166      ! 
    166167   END SUBROUTINE iom_swap 
    167168 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90

    r3680 r4220  
    4444 
    4545   INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100   !: maximum number of simultaneously opened file 
    46    INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 360  !: maximum number of variables in one file 
     46   INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 600  !: maximum number of variables in one file 
    4747   INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable 
    4848   INTEGER, PARAMETER, PUBLIC ::   jpmax_digits =  5   !: maximum number of digits for the cpu number in the file name 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r3680 r4220  
    2525   USE domvvl          ! variable volume 
    2626   USE divcur          ! hor. divergence and curl      (div & cur routines) 
     27   USE sbc_ice, ONLY : lk_lim3 
    2728 
    2829   IMPLICIT NONE 
     
    132133                     CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       ) 
    133134#endif 
     135                  IF( lk_lim3 ) THEN 
     136                     CALL iom_rstput( kt, nitrst, numrow, 'iatte'  , iatte     ) !clem modif 
     137                     CALL iom_rstput( kt, nitrst, numrow, 'oatte'  , oatte     ) !clem modif 
     138                  ENDIF 
    134139      IF( kt == nitrst ) THEN 
    135140         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
     
    232237      ENDIF 
    233238      ! 
     239      IF( lk_lim3 ) THEN  
     240         CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 
     241         CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif 
     242      ENDIF 
     243      ! 
    234244   END SUBROUTINE rst_read 
    235245 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r4155 r4220  
    4040   LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS       
    4141   LOGICAL , PUBLIC ::   ln_apr_dyn  = .FALSE.   !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
    42    LOGICAL , PUBLIC ::   ln_icebergs = .FALSE.   !: Icebergs 
    4342   INTEGER , PUBLIC ::   nn_ice      = 0         !: flag for ice in the surface boundary condition (=0/1/2/3) 
    4443   INTEGER , PUBLIC ::   nn_ice_embd = 0         !: flag for levitating/embedding sea-ice in the ocean 
     
    5352   LOGICAL , PUBLIC ::   ln_cdgw     = .FALSE.   !: true if neutral drag coefficient from wave model 
    5453   LOGICAL , PUBLIC ::   ln_sdw      = .FALSE.   !: true if 3d stokes drift from wave model 
     54   ! 
     55   LOGICAL , PUBLIC ::   ln_icebergs = .FALSE.   !: Icebergs 
    5556   ! 
    5657   CHARACTER (len=8), PUBLIC :: cn_iceflx = 'none' !: Flux handling over ice categories 
     
    8182   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s] jpi,jpj,jpk 
    8283   !! 
    83    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   oatte, iatte      !: clem attenuation coef of the input solar flux [unitless] 
    8484   !! 
    8585   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip           !: total precipitation                          [Kg/m2/s] 
     
    126126         ! 
    127127      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
    128          &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,     & 
    129          &      iatte(jpi,jpj) , oatte    (jpi,jpj)                              , STAT=ierr(3) ) 
     128         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
    130129         ! 
    131130      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r4099 r4220  
    388388      ! 
    389389      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used 
    390          srcv(jpr_itz1:jpr_itz2)%laction = .FALSE.    ! ice components not received (itx1 and ity1 used later) 
     390         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received 
    391391         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation 
    392392         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp. 
     
    407407      SELECT CASE( TRIM( sn_rcv_emp%cldes ) ) 
    408408      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE.  
    409       CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 
     409      CASE( 'conservative'  ) 
     410         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 
     411         IF ( k_ice <= 1 )  srcv(jpr_ivep)%laction = .FALSE. 
    410412      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE. 
    411413      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) 
     
    465467         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' ) 
    466468      !                                                      ! ------------------------- ! 
    467       !                                                      !    Ice Qsr penetration    !    
    468       !                                                      ! ------------------------- ! 
    469       ! fraction of net shortwave radiation which is not absorbed in the thin surface layer  
    470       ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
    471       ! Coupled case: since cloud cover is not received from atmosphere  
    472       !               ===> defined as constant value -> definition done in sbc_cpl_init 
    473       IF ( ALLOCATED (fr1_i0)) fr1_i0 (:,:) = 0.18 
    474       IF ( ALLOCATED (fr2_i0)) fr2_i0 (:,:) = 0.82 
    475       !                                                      ! ------------------------- ! 
    476469      !                                                      !      10m wind module      !    
    477470      !                                                      ! ------------------------- ! 
     
    508501      ! Allocate taum part of frcv which is used even when not received as coupling field 
    509502      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) ) 
     503      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE. 
     504      IF( k_ice /= 0 ) THEN 
     505         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jn)%nct) ) 
     506         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jn)%nct) ) 
     507      END IF 
    510508 
    511509      ! ================================ ! 
     
    13311329      END SELECT 
    13321330 
     1331      !    Ice Qsr penetration used (only?)in lim2 or lim3  
     1332      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer  
     1333      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
     1334      ! Coupled case: since cloud cover is not received from atmosphere  
     1335      !               ===> defined as constant value -> definition done in sbc_cpl_init 
     1336      fr1_i0(:,:) = 0.18 
     1337      fr2_i0(:,:) = 0.82 
     1338 
     1339 
    13331340      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr ) 
    13341341      ! 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4099 r4220  
    153153      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    154154 
    155       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    156          CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
    157       END IF 
    158       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    159          CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    160       ENDIF 
    161  
     155#if defined key_coupled 
     156      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
     157      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
     158         &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     159#endif 
    162160 
    163161      IF( kt == nit000 ) THEN 
     
    169167         ! 
    170168         IF( ln_nicep ) THEN      ! control print at a given point 
    171             jiindx = 15   ;   jjindx = 46 
     169            jiindx = 6   ;   jjindx = 47 
    172170            WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    173171         ENDIF 
     
    193191         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
    194192          
     193#if defined key_coupled 
    195194         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    196195            ! 
     
    200199            ! 
    201200         ENDIF 
    202                                                      ! Bulk formulea - provides the following fields: 
     201#endif 
     202                                               ! Bulk formulea - provides the following fields: 
    203203         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
    204204         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2] 
     
    238238 
    239239         ! Average over all categories 
     240#if defined key_coupled 
    240241         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    241242 
     
    267268            END IF 
    268269         END IF 
    269  
     270#endif 
    270271         !                                           !----------------------! 
    271272         !                                           ! LIM-3  time-stepping ! 
     
    283284         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content 
    284285         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content 
    285  
     286         ! 
     287         old_u_ice(:,:) = u_ice(:,:) 
     288         old_v_ice(:,:) = v_ice(:,:) 
    286289         !                                           ! intialisation to zero    !!gm is it truly necessary ??? 
    287290         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
     
    292295         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
    293296         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
     297         ! 
     298         d_u_ice_dyn(:,:) = 0._wp 
     299         d_v_ice_dyn(:,:) = 0._wp 
    294300         ! 
    295301         sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp 
     
    358364                          CALL lim_update2                 ! Global variables update 
    359365#if defined key_bdy 
     366                          CALL lim_var_glo2eqv            ! 
    360367                          CALL bdy_ice_lim( kt )          ! clem modif: bdy ice 
    361368#endif 
     
    373380                          CALL lim_wri( 1  )              ! Ice outputs  
    374381!clem # endif 
    375          IF( kt == nit000 )   CALL iom_close( numrir )  ! clem: close input ice restart file 
     382         IF( kt == nit000 .AND. ln_rstart )   & 
     383            &             CALL iom_close( numrir )        ! clem: close input ice restart file 
     384         ! 
    376385         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file  
    377386                          CALL lim_var_glo2eqv            ! ??? 
     
    392401      ! 
    393402      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    394       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    395          CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
    396       END IF 
    397       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    398          CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    399       ENDIF 
     403 
     404#if defined key_coupled 
     405      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
     406      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
     407         &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     408#endif 
    400409      ! 
    401410      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4155 r4220  
    106106          nn_ice      =   0 
    107107      ENDIF 
    108        
     108      
    109109      IF(lwp) THEN               ! Control print 
    110110         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)' 
     
    129129      ENDIF 
    130130 
    131       !   Flux handling over ice categories  
     131      !   Flux handling over ice categories 
     132#if defined key_coupled  
    132133      SELECT CASE ( TRIM (cn_iceflx)) 
    133134      CASE ('ave') 
     
    143144      IF(lwp) WRITE(numout,*) '              Fluxes averaged over all ice categories         ln_iceflx_ave    = ', ln_iceflx_ave 
    144145      IF(lwp) WRITE(numout,*) '              Fluxes distributed linearly over ice categories ln_iceflx_linear = ', ln_iceflx_linear 
     146#endif 
    145147      ! 
    146148      !                              ! allocate sbc arrays 
     
    182184      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   & 
    183185         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' ) 
    184  
     186#if defined key_coupled 
    185187      IF( ln_iceflx_ave .AND. ln_iceflx_linear ) & 
    186188         &   CALL ctl_stop( ' ln_iceflx_ave and ln_iceflx_linear options are not compatible' ) 
    187  
    188189      IF( ( nn_ice ==3 .AND. lk_cpl) .AND. .NOT. ( ln_iceflx_ave .OR. ln_iceflx_linear ) ) & 
    189190         &   CALL ctl_stop( ' With lim3 coupled, either ln_iceflx_ave or ln_iceflx_linear must be set to .TRUE.' ) 
    190        
     191#endif       
    191192      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag 
    192193 
     
    320321      CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model 
    321322      !is it useful? 
    322       !CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model 
     323      CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model 
    323324      END SELECT                                               
    324325 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4038 r4220  
    3232   USE wrk_nemo       ! Memory Allocation 
    3333   USE timing         ! Timing 
     34   USE sbc_ice, ONLY : lk_lim3 
    3435 
    3536 
     
    163164         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution 
    164165         ! clem: store attenuation coefficient of the first ocean level 
    165          IF ( ln_qsr_ice ) THEN 
     166         IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    166167            DO jj = 1, jpj 
    167168               DO ji = 1, jpi 
     
    232233               END DO 
    233234               ! clem: store attenuation coefficient of the first ocean level 
    234                IF ( ln_qsr_ice ) THEN 
     235               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    235236                  DO jj = 1, jpj 
    236237                     DO ji = 1, jpi 
     
    256257               END DO 
    257258               ! clem: store attenuation coefficient of the first ocean level 
    258                IF ( ln_qsr_ice ) THEN 
     259               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    259260                  oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    260261                  iatte(:,:) = oatte(:,:) 
     
    280281               END DO 
    281282               ! clem: store attenuation coefficient of the first ocean level 
    282                IF ( ln_qsr_ice ) THEN 
     283               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    283284                  DO jj = 1, jpj 
    284285                     DO ji = 1, jpi 
     
    299300               END DO 
    300301               ! clem: store attenuation coefficient of the first ocean level 
    301                IF ( ln_qsr_ice ) THEN 
     302               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
    302303                  oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    303304                  iatte(:,:) = oatte(:,:) 
     
    319320         ! 
    320321      ENDIF 
    321       ! clem: store attenuation coefficient of the first ocean level 
    322       !IF (ln_traqsr) THEN 
    323       !   DO jj = 1, jpj 
    324       !      DO ji = 1, jpi 
    325       !         IF ( qsr(ji,jj) /= 0._wp ) THEN 
    326       !            oatte(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 
    327       !            iatte(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 
    328       !         ENDIF 
    329       !      END DO 
    330       !   END DO 
    331       !END IF 
    332322      ! 
    333323      IF( lrst_oce ) THEN   !                  Write in the ocean restart file 
     
    375365      !!---------------------------------------------------------------------- 
    376366      ! 
    377       INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     367      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    378368      INTEGER  ::   irgb, ierror, ioptio, nqsr   ! local integer 
    379369      REAL(wp) ::   zz0, zc0  , zc1, zcoef       ! local scalars 
     
    393383      ! 
    394384      ! clem init for oatte and iatte 
    395       oatte(:,:) = 1._wp 
    396       iatte(:,:) = 1._wp 
     385      IF( .NOT. ln_rstart ) THEN 
     386         oatte(:,:) = 1._wp 
     387         iatte(:,:) = 1._wp 
     388      ENDIF 
    397389      ! 
    398390      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
     
    423415         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
    424416         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1 
     417         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice     
    425418      ENDIF 
    426419 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r3625 r4220  
    5555   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_fmass       !: time evolution of mass of snow+ice               [Kg/m2/s] 
    5656 
     57   !! arrays related to penetration of solar fluxes to calculate the heat budget for sea ice 
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   oatte, iatte       !: attenuation coef of the input solar flux [unitless] 
     59 
    5760   !!---------------------------------------------------------------------- 
    5861   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     
    6669      !!                   ***  FUNCTION oce_alloc  *** 
    6770      !!---------------------------------------------------------------------- 
    68       INTEGER :: ierr(3) 
     71      INTEGER :: ierr(4) 
    6972      !!---------------------------------------------------------------------- 
    7073      ! 
     
    8790         &      gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) ) 
    8891         ! 
    89       ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             & 
    90          &      snwice_fmass(jpi,jpj), STAT= ierr(3) ) 
     92      ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) ) 
     93         ! 
     94      ALLOCATE( iatte(jpi,jpj) , oatte(jpi,jpj) , STAT=ierr(4) ) 
    9195         ! 
    9296      oce_alloc = MAXVAL( ierr ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4155 r4220  
    9090      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9191                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     92 
    9293      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp ) 
    9394      IF( lk_tide    )   CALL sbc_tide( kstp ) 
     
    108109      ! 
    109110      !  VERTICAL PHYSICS 
    110       ! bg jchanut tschanges 
    111       ! One need bottom friction parameter in ssh_wzv routine with time splitting. 
    112       ! The idea could be to move the call below before ssh_wzv. However, "now" scale factors 
    113       ! at U-V points (which are set thanks to sshu_n, sshv_n) are actually available in sshwzv. 
    114       ! These are needed for log bottom friction... 
    115 #if ! defined key_dynspg_ts 
    116111                         CALL zdf_bfr( kstp )         ! bottom friction 
    117 #endif 
    118       ! end jchanut tschanges 
    119112 
    120113      !                                               ! Vertical eddy viscosity and diffusivity coefficients 
     
    214207            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    215208 
    216       ELSE    
    217                                                ! centered hpg  (eos then time stepping) 
    218       ! bg jchanut tschanges 
    219 #if ! defined key_dynspg_ts 
    220       ! eos already called 
     209      ELSE                                                  ! centered hpg  (eos then time stepping) 
    221210                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation 
    222211         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative 
    223212            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    224 #endif 
    225       ! end jchanut tschanges 
    226213         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    227214                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
     
    231218      ! Dynamics                                    (tsa used as workspace) 
    232219      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    233       ! bg jchanut tschanges 
    234 #if defined key_dynspg_ts       
    235 ! revert to previously computed tendencies: 
    236 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 
    237                                ua(:,:,:) = ua_bak(:,:,:)             
    238                                va(:,:,:) = va_bak(:,:,:) 
    239                                CALL dyn_bfr( kstp )         ! bottom friction 
    240                                CALL dyn_zdf( kstp )         ! vertical diffusion 
    241 #else 
    242       ! end jchanut tschanges 
    243220                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    244221                               va(:,:,:) = 0.e0 
     
    260237                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    261238                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    262       ! bg jchanut tschanges 
    263 #endif 
    264       ! end jchanut tschanges 
    265239                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
    266240 
Note: See TracChangeset for help on using the changeset viewer.