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 4649 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 – NEMO

Ignore:
Timestamp:
2014-05-27T11:28:12+02:00 (10 years ago)
Author:
clem
Message:

finalizing LIM3 heat budget conservation + multiple minor bugs corrections

File:
1 edited

Legend:

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

    r4634 r4649  
    3232   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3333 
    34    REAL(wp) ::   epsi10      = 1.e-10_wp    ! 
     34   REAL(wp) ::   epsi10 = 1.e-10_wp    ! 
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     
    9292      !!           (04-2007) Energy conservation tested by M. Vancoppenolle 
    9393      !!------------------------------------------------------------------ 
    94       INTEGER , INTENT (in) ::   kideb   ! Start point on which the  the computation is applied 
    95       INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied 
    96       INTEGER , INTENT (in) ::   jl      ! Category number 
     94      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
     95      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    9796 
    9897      !! * Local variables 
     
    103102      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    104103      INTEGER ::   minnumeqmin, maxnumeqmax 
    105       INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation 
    106       INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    107       INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     104      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
     105      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
     106      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    108107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    109108      REAL(wp) ::   zg1       =  2._wp        ! 
     
    115114      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    116115      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    117       REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point 
    118       REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure ) 
    119       REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration 
    120       REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness 
    121       REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness 
    122       REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface 
    123       REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function 
    124       REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function 
    125       REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature 
    126       REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4) 
    127       REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice 
    128       REAL(wp), DIMENSION(kiut) ::   zihic, zhsu 
    129       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
    130       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
    131       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
    132       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
    133       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice 
    134       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
    135       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    136       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat 
    137       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice 
    138       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
    139       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
    140       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
    141       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow 
    142       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence 
    143       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow 
    144       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow 
    145       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term 
    146       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term 
    147       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis 
    148       REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms 
    149       REAL(wp) ::   ztemp   ! local scalar  
     116      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point 
     117      REAL(wp), POINTER, DIMENSION(:) ::   ztsuold     ! old surface temperature (before the iterative procedure ) 
     118      REAL(wp), POINTER, DIMENSION(:) ::   ztsuoldit   ! surface temperature at previous iteration 
     119      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
     120      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
     121      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface 
     122      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function 
     123      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function 
     124      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature 
     125      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4) 
     126      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice 
     127      REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu 
     128      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity 
     129      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice 
     130      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice 
     131      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   ztiold      ! Old temperature in the ice 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice 
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     135      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow 
     140      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s       ! Eta factor in the snow 
     141      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp      ! Temporary temperature in the snow to check the convergence 
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsold       ! Temporary temperature in the snow 
     143      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s          ! Vertical cotes of the layers in the snow 
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm   ! Independent term 
     145      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis   ! temporary independent term 
     146      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
     147      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms 
    150148      !!------------------------------------------------------------------      
    151149      !  
     150      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
     151      CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 
     152      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
     153      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
     154      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0) 
     155      CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  ) 
     156      CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) 
     157 
    152158      !------------------------------------------------------------------------------! 
    153159      ! 1) Initialization                                                            ! 
     
    318324            DO layer = 1, nlay_i-1 
    319325               DO ji = kideb , kiut 
    320                   ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt 
    321                   ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    322                      &                                   / MIN( -2.0_wp * epsi10, ztemp )   & 
    323                      &                       - 0.0055_wp * ztemp 
     326                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
     327                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
     328                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    324329                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
    325330               END DO 
    326331            END DO 
    327332            DO ji = kideb , kiut 
    328                ztemp = t_bo_b(ji) - rtt 
    329                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp )   & 
    330                   &                        - 0.011_wp * ztemp   
     333               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
     334                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
    331335               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    332336            END DO 
     
    395399         ! 
    396400         DO ji = kideb , kiut 
    397  
    398401            ! update of the non solar flux according to the update in T_su 
    399402            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) 
     
    403406               + qns_ice_1d(ji)                  ! non solar total flux  
    404407            ! (LWup, LWdw, SH, LH) 
    405  
    406             ! heat flux used to change surface temperature 
    407             !hfx_tot_1d(ji) = hfx_tot_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) * a_i_b(ji) 
    408408         END DO 
    409409 
     
    721721      DO ji = kideb, kiut 
    722722         IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC 
    723             hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
     723            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    724724         ELSE                                    ! case T_su = 0degC 
    725             hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
     725            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    726726         ENDIF 
    727727      END DO 
    728728 
    729729      ! 
     730      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
     731      CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 
     732      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
     733      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
     734      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 
     735      CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 
     736      CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 
     737 
    730738   END SUBROUTINE lim_thd_dif 
    731739 
Note: See TracChangeset for help on using the changeset viewer.