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 14009 for NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/icethd_zdf_bl99.F90 – NEMO

Ignore:
Timestamp:
2020-12-02T15:42:07+01:00 (3 years ago)
Author:
gsamson
Message:

dev_r2052_ENHANCE-09_rbourdal_massfluxconvection update with trunk@r14008

Location:
NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette_MPI3_LoopFusion@13943         sette 
         10^/utils/CI/sette_wave@13990         sette 
  • NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/icethd_zdf_bl99.F90

    r13472 r14009  
    109109      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function 
    110110      ! 
    111       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    112       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    113       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    114       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
    115       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    116       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    117       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    119       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    120       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    121       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    122       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    123       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    124       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    125       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    126       REAL(wp), DIMENSION(jpij)            ::   zkappa_comb ! Combined snow and ice surface conductivity 
    127       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    128       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    129       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    130       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    131       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    132       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    133       REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow  
    134       REAL(wp), DIMENSION(jpij)            ::   isnow       ! snow presence (1) or not (0)  
    135       REAL(wp), DIMENSION(jpij)            ::   isnow_comb  ! snow presence for met-office  
     111      REAL(wp), DIMENSION(jpij       )   ::   ztsuold     ! Old surface temperature in the ice 
     112      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztiold      ! Old temperature in the ice 
     113      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsold      ! Old temperature in the snow 
     114      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Temporary temperature in the ice to check the convergence 
     115      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow to check the convergence 
     116      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
     117      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i_cp ! copy 
     118      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
     119      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
     120      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
     121      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
     122      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
     123      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
     124      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
     125      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zeta_s      ! Eta factor in the snow 
     126      REAL(wp), DIMENSION(jpij)          ::   zkappa_comb ! Combined snow and ice surface conductivity 
     127      REAL(wp), DIMENSION(jpij)          ::   zq_ini      ! diag errors on heat 
     128      REAL(wp), DIMENSION(jpij)          ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     129      REAL(wp), DIMENSION(jpij)          ::   za_s_fra    ! ice fraction covered by snow  
     130      REAL(wp), DIMENSION(jpij)          ::   isnow       ! snow presence (1) or not (0)  
     131      REAL(wp), DIMENSION(jpij)          ::   isnow_comb  ! snow presence for met-office  
     132      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindterm    ! 'Ind'ependent term 
     133      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindtbis    ! Temporary 'ind'ependent term 
     134      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     135      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) ::   ztrid       ! Tridiagonal system terms 
    136136      ! 
    137137      ! Mono-category 
     
    533533            ! Solve the tridiagonal system with Gauss elimination method. 
    534534            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    535             jm_maxt = 0 
    536             jm_mint = nlay_i+5 
    537             DO ji = 1, npti 
    538                jm_mint = MIN(jm_min(ji),jm_mint) 
    539                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    540             END DO 
    541  
    542             DO jk = jm_mint+1, jm_maxt 
    543                DO ji = 1, npti 
    544                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     535!!$            jm_maxt = 0 
     536!!$            jm_mint = nlay_i+5 
     537!!$            DO ji = 1, npti 
     538!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     539!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     540!!$            END DO 
     541!!$            !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 
     542!!$             
     543!!$            DO jk = jm_mint+1, jm_maxt 
     544!!$               DO ji = 1, npti 
     545!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     546!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     547!!$                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     548!!$               END DO 
     549!!$            END DO 
     550            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     551            DO ji = 1, npti 
     552               jm_mint = jm_min(ji) 
     553               jm_maxt = jm_max(ji) 
     554               DO jm = jm_mint+1, jm_maxt 
    545555                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    546556                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     
    564574            END DO 
    565575 
     576            ! snow temperatures       
    566577            DO ji = 1, npti 
    567578               ! Variables used after iterations 
    568579               ! Value must be frozen after convergence for MPP independance reason 
    569                IF ( .NOT. l_T_converged(ji) ) THEN 
    570                   ! snow temperatures       
    571                   IF( h_s_1d(ji) > 0._wp ) THEN 
    572                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    573                   ENDIF 
    574                   ! surface temperature 
     580               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     581                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     582            END DO 
     583            !!clem SNWLAY 
     584            DO jm = nlay_s, 2, -1 
     585               DO ji = 1, npti 
     586                  jk = jm - 1 
     587                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     588                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     589               END DO 
     590            END DO 
     591             
     592            ! surface temperature 
     593            DO ji = 1, npti 
     594               IF( .NOT. l_T_converged(ji) ) THEN 
    575595                  ztsub(ji) = t_su_1d(ji) 
    576596                  IF( t_su_1d(ji) < rt0 ) THEN 
    577                      t_su_1d(ji) = (  zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    578                         &           ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     597                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     598                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    579599                  ENDIF 
    580600               ENDIF 
    581601            END DO 
    582             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
    583602            ! 
    584603            !-------------------------------------------------------------- 
     
    727746            ! Solve the tridiagonal system with Gauss elimination method. 
    728747            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    729             jm_maxt = 0 
    730             jm_mint = nlay_i+5 
    731             DO ji = 1, npti 
    732                jm_mint = MIN(jm_min(ji),jm_mint) 
    733                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    734             END DO 
    735              
    736             DO jk = jm_mint+1, jm_maxt 
    737                DO ji = 1, npti 
    738                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     748!!$            jm_maxt = 0 
     749!!$            jm_mint = nlay_i+5 
     750!!$            DO ji = 1, npti 
     751!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     752!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     753!!$            END DO 
     754!!$             
     755!!$            DO jk = jm_mint+1, jm_maxt 
     756!!$               DO ji = 1, npti 
     757!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     758!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     759!!$                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1) 
     760!!$               END DO 
     761!!$            END DO 
     762            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     763            DO ji = 1, npti 
     764               jm_mint = jm_min(ji) 
     765               jm_maxt = jm_max(ji) 
     766               DO jm = jm_mint+1, jm_maxt 
    739767                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    740                   zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)  / zdiagbis(ji,jm-1) 
     768                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
    741769               END DO 
    742770            END DO 
    743              
     771 
    744772            ! ice temperatures 
    745773            DO ji = 1, npti 
     
    761789            ! snow temperatures       
    762790            DO ji = 1, npti 
    763                ! Variable used after iterations 
     791               ! Variables used after iterations 
    764792               ! Value must be frozen after convergence for MPP independance reason 
    765                IF ( .NOT. l_T_converged(ji) ) THEN 
    766                   IF( h_s_1d(ji) > 0._wp ) THEN 
    767                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    768                   ENDIF 
    769                ENDIF 
    770             END DO 
    771             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
     793               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     794                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     795            END DO 
     796            !!clem SNWLAY 
     797            DO jm = nlay_s, 2, -1 
     798               DO ji = 1, npti 
     799                  jk = jm - 1 
     800                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     801                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     802               END DO 
     803            END DO 
    772804            ! 
    773805            !-------------------------------------------------------------- 
     
    923955         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    924956         IF( h_s_1d(ji) >= zhs_ssl ) THEN 
    925             t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1)   & 
    926                &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 
     957            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s)   & 
     958               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1)      ) & 
    927959               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i & 
    928960               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 
Note: See TracChangeset for help on using the changeset viewer.