Ignore:
Timestamp:
2015-02-02T11:28:50+01:00 (6 years ago)
Author:
vancop
Message:

new itd boundaries, namelist changes, mono-category and comments

File:
1 edited

Legend:

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

    r5047 r5048  
    100100      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    101101      INTEGER ::   minnumeqmin, maxnumeqmax 
     102       
    102103      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
    103104      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
    104105      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     106       
    105107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    106108      REAL(wp) ::   zg1       =  2._wp        ! 
     
    113115      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    114116      REAL(wp) ::   zhsu 
    115       REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point 
    116       REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    117       REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration 
    118       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
    119       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    120       REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface 
    121       REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function 
    122       REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function 
    123       REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature 
    124       REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4) 
    125       REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice 
    126       REAL(wp), POINTER, DIMENSION(:) ::   zihic 
    127       REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity 
    128       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice 
    130       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow 
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zswiterm    ! Independent term 
    144       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitbis    ! temporary independent term 
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
    146       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms 
     117       
     118      REAL(wp), POINTER, DIMENSION(:)     ::   ztfs        ! ice melting point 
     119      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
     120      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration 
     121      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness 
     122      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness 
     123      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface 
     124      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function 
     125      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function 
     126      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature 
     127      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4) 
     128      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice 
     129      REAL(wp), POINTER, DIMENSION(:)     ::   zihic 
     130       
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice 
     133      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice 
     134      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice 
     135      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice 
     136      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice 
     137      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat 
     139      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice 
     140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow 
     141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow 
     142      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow 
     143      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow 
     144      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
     145      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow 
     146      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow 
     147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term 
     148      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term 
     149      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms 
     151       
    147152      ! diag errors on heat 
    148       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 
     153      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err 
     154       
     155      ! Mono-category 
     156      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done 
     157      REAL(wp)                            :: zratio_s      ! dummy factor 
     158      REAL(wp)                            :: zratio_i      ! dummy factor 
     159      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation 
     160      REAL(wp)                            :: zhe           ! dummy factor 
     161      REAL(wp)                            :: zswitch       ! dummy switch 
     162      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity 
     163      REAL(wp)                            :: zfac          ! dummy factor 
     164      REAL(wp)                            :: zihe          ! dummy factor 
     165      REAL(wp)                            :: zheshth       ! dummy factor 
     166       
     167      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat 
     168       
    149169      !!------------------------------------------------------------------      
    150170      !  
    151171      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
    152172      CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    153       CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic ) 
     173      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    154174      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
    155175      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 
    156       CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis  ) 
     176      CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  ) 
    157177      CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 
    158178 
     
    200220      ! 
    201221      !------------------------------------------------------------------------------| 
    202       ! 2) Radiations                                                                | 
     222      ! 2) Radiation                                              | 
    203223      !------------------------------------------------------------------------------| 
    204224      ! 
     
    213233      ! zftrice = io.qsr_ice       is below the surface  
    214234      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
     235      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    215236      zhsu = 0.1_wp ! threshold for the computation of i0 
    216          !fr1_i0_1d = i0 for a thin ice surface 
    217          !fr1_i0_2d = i0 for a thick ice surface 
    218237      DO ji = kideb , kiut 
    219238         ! switches 
     
    339358            END DO 
    340359         ENDIF 
    341          ! 
    342          !------------------------------------------------------------------------------| 
    343          !  5) kappa factors                                                            | 
    344          !------------------------------------------------------------------------------| 
    345          ! 
     360          
     361         ! 
     362         !------------------------------------------------------------------------------| 
     363         !  6) G(he) - enhancement of thermal conductivity in mono-category case        | 
     364         !------------------------------------------------------------------------------| 
     365         ! 
     366         ! Computation of effective thermal conductivity G(h) 
     367         ! Used in mono-category case only to simulate an ITD implicitly 
     368         ! Fichefet and Morales Maqueda, JGR 1997 
     369 
     370         zghe(:) = 0._wp 
     371 
     372         SELECT CASE ( nn_monocat ) 
     373 
     374         CASE (0,2,4) 
     375 
     376            zghe(kideb:kiut) = 1._wp 
     377 
     378         CASE (1,3) ! LIM3 
     379 
     380            zepsilon = 0.1 
     381            zh_thres = EXP( 1._wp ) * zepsilon / 2. 
     382 
     383            DO ji = kideb, kiut 
     384    
     385               ! Mean sea ice thermal conductivity 
     386               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp) 
     387 
     388               ! Effective thickness he (zhe) 
     389               zfac     = 1._wp / ( rcdsn + zkimean ) 
     390               zratio_s = rcdsn   * zfac 
     391               zratio_i = zkimean * zfac 
     392               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     393 
     394               ! G(he) 
     395               zswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
     396               zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) ) 
     397    
     398               ! Impose G(he) < 2. 
     399               zghe(ji) = MIN( zghe(ji), 2.0 ) 
     400 
     401            END DO 
     402 
     403         END SELECT 
     404 
     405         ! 
     406         !------------------------------------------------------------------------------| 
     407         !  7) kappa factors                                                            | 
     408         !------------------------------------------------------------------------------| 
     409         ! 
     410         !--- Snow 
    346411         DO ji = kideb, kiut 
    347  
    348             !-- Snow kappa factors 
    349             zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
    350             zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
     412            zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
     413            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac 
     414            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac 
    351415         END DO 
    352416 
    353417         DO jk = 1, nlay_s-1 
    354418            DO ji = kideb , kiut 
    355                zkappa_s(ji,jk)  = 2.0 * rcdsn / & 
    356                   MAX(epsi10,2.0*zh_s(ji)) 
    357             END DO 
    358          END DO 
    359  
     419               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) ) 
     420            END DO 
     421         END DO 
     422 
     423         !--- Ice 
    360424         DO jk = 1, nlay_i-1 
    361425            DO ji = kideb , kiut 
    362                !-- Ice kappa factors 
    363                zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ & 
    364                   MAX(epsi10,2.0*zh_i(ji))  
    365             END DO 
    366          END DO 
    367  
    368          DO ji = kideb , kiut 
    369             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
    370             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    371             !-- Interface 
    372             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    373                (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    374             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
    375                + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 
    376          END DO 
    377          ! 
    378          !------------------------------------------------------------------------------| 
    379          ! 6) Sea ice specific heat, eta factors                                        | 
     426               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) ) 
     427            END DO 
     428         END DO 
     429 
     430         !--- Snow-ice interface 
     431         DO ji = kideb , kiut 
     432            zfac                  = 1./ MAX( epsi10 , zh_i(ji) ) 
     433            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
     434            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
     435            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
     436           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
     437            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp ) 
     438         END DO 
     439 
     440         ! 
     441         !------------------------------------------------------------------------------| 
     442         ! 8) Sea ice specific heat, eta factors                                        | 
    380443         !------------------------------------------------------------------------------| 
    381444         ! 
     
    383446            DO ji = kideb , kiut 
    384447               ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    385                zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 
    386                   MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 
    387                zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 
    388                   epsi10) 
     448               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 ) 
     449               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 ) 
    389450            END DO 
    390451         END DO 
     
    396457            END DO 
    397458         END DO 
    398          ! 
    399          !------------------------------------------------------------------------------| 
    400          ! 7) surface flux computation                                                  | 
     459 
     460         ! 
     461         !------------------------------------------------------------------------------| 
     462         ! 9) surface flux computation                                                  | 
    401463         !------------------------------------------------------------------------------| 
    402464         ! 
     
    418480         ! 
    419481         !------------------------------------------------------------------------------| 
    420          ! 8) tridiagonal system terms                                                  | 
     482         ! 10) tridiagonal system terms                                                 | 
    421483         !------------------------------------------------------------------------------| 
    422484         ! 
     
    433495               ztrid(ji,numeq,2) = 0. 
    434496               ztrid(ji,numeq,3) = 0. 
    435                zswiterm(ji,numeq)= 0. 
    436                zswitbis(ji,numeq)= 0. 
     497               zindterm(ji,numeq)= 0. 
     498               zindtbis(ji,numeq)= 0. 
    437499               zdiagbis(ji,numeq)= 0. 
    438500            ENDDO 
     
    446508                  zkappa_i(ji,jk)) 
    447509               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk) 
    448                zswiterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
     510               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
    449511                  zradab_i(ji,jk) 
    450512            END DO 
     
    458520               +  zkappa_i(ji,nlay_i-1) ) 
    459521            ztrid(ji,numeq,3)  =  0.0 
    460             zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
     522            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
    461523               ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 
    462524               *  t_bo_1d(ji) )  
     
    478540                     zkappa_s(ji,jk) ) 
    479541                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    480                   zswiterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
     542                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
    481543                     zradab_s(ji,jk) 
    482544               END DO 
     
    485547               IF ( nlay_i.eq.1 ) THEN 
    486548                  ztrid(ji,nlay_s+2,3)    =  0.0 
    487                   zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
     549                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
    488550                     t_bo_1d(ji)  
    489551               ENDIF 
     
    502564                  ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 
    503565                  ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 
    504                   zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji) 
     566                  zindterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji) 
    505567 
    506568                  !!first layer of snow equation 
     
    508570                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 
    509571                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    510                   zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
     572                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
    511573 
    512574               ELSE  
     
    525587                     zkappa_s(ji,0) * zg1s ) 
    526588                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    527                   zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
     589                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
    528590                     ( zradab_s(ji,1) +                         & 
    529591                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     
    549611                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
    550612                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    551                   zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
     613                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
    552614 
    553615                  !!first layer of ice equation 
     
    556618                     + zkappa_i(ji,0) * zg1 ) 
    557619                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)   
    558                   zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
     620                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
    559621 
    560622                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    569631                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    570632 
    571                      zswiterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
     633                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
    572634                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 
    573635                  ENDIF 
     
    589651                     zg1)   
    590652                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    591                   zswiterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
     653                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
    592654                     zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    593655 
     
    598660                        zkappa_i(ji,1)) 
    599661                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    600                      zswiterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
     662                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
    601663                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 
    602664                        + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 
     
    610672         ! 
    611673         !------------------------------------------------------------------------------| 
    612          ! 9) tridiagonal system solving                                                | 
     674         ! 11) tridiagonal system solving                                               | 
    613675         !------------------------------------------------------------------------------| 
    614676         ! 
     
    622684 
    623685         DO ji = kideb , kiut 
    624             zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji)) 
     686            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
    625687            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
    626688            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin) 
     
    633695               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 
    634696                  ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 
    635                zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 
    636                   zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
     697               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 
     698                  zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
    637699            END DO 
    638700         END DO 
     
    640702         DO ji = kideb , kiut 
    641703            ! ice temperatures 
    642             t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
     704            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
    643705         END DO 
    644706 
     
    646708            DO ji = kideb , kiut 
    647709               jk    =  numeq - nlay_s - 1 
    648                t_i_1d(ji,jk)  =  (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 
     710               t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 
    649711                  t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 
    650712            END DO 
     
    654716            ! snow temperatures       
    655717            IF (ht_s_1d(ji).GT.0._wp) & 
    656                t_s_1d(ji,nlay_s)     =  (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
     718               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    657719               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 
    658720               *        MAX(0.0,SIGN(1.0,ht_s_1d(ji)))  
     
    662724            ztsubit(ji) = t_su_1d(ji) 
    663725            IF( t_su_1d(ji) < ztfs(ji) ) & 
    664                t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
     726               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
    665727               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    666728         END DO 
    667729         ! 
    668730         !-------------------------------------------------------------------------- 
    669          !  10) Has the scheme converged ?, end of the iterative procedure         | 
     731         !  12) Has the scheme converged ?, end of the iterative procedure         | 
    670732         !-------------------------------------------------------------------------- 
    671733         ! 
     
    709771      ! 
    710772      !-------------------------------------------------------------------------! 
    711       !   11) Fluxes at the interfaces                                          ! 
     773      !   13) Fluxes at the interfaces                                          ! 
    712774      !-------------------------------------------------------------------------! 
    713775      DO ji = kideb, kiut 
     
    771833      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    772834      CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    773       CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic ) 
     835      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    774836      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 
    775837         &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
    776838      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
    777       CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 
     839      CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 
    778840      CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 
    779841      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 
Note: See TracChangeset for help on using the changeset viewer.