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 10425 for NEMO/trunk/src/ICE/icethd_zdf_bl99.F90 – NEMO

Ignore:
Timestamp:
2018-12-19T22:54:16+01:00 (5 years ago)
Author:
smasson
Message:

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icethd_zdf_bl99.F90

    r10069 r10425  
    8383      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation 
    8484      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation 
    85       ! 
     85 
     86      LOGICAL, DIMENSION(jpij) ::   l_T_converged   ! true when T converges (per grid point) 
     87! 
    8688      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    8789      REAL(wp) ::   zg1       =  2._wp        ! 
     
    113115      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    114116      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
     117      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy 
    115118      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    116119      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
     
    201204      ! 
    202205      iconv    = 0          ! number of iterations 
    203       zdti_max = 1000._wp   ! maximal value of error on all points 
    204       ! 
     206      ! 
     207      l_T_converged(:) = .FALSE. 
    205208      !                                                          !============================! 
    206       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
     209      ! Convergence calculated until all sub-domain grid points have converged 
     210      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 
     211      ! but values are not taken into account (results independant of MPI partitioning) 
     212      ! 
     213      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    207214         !                                                       !============================! 
    208215         iconv = iconv + 1 
     
    217224            ! 
    218225            DO ji = 1, npti 
    219                ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    220                ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
     226               ztcond_i_cp(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     227               ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    221228            END DO 
    222229            DO jk = 1, nlay_i-1 
    223230               DO ji = 1, npti 
    224                   ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    225                      &                       MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     231                  ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
     232                     &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    226233               END DO 
    227234            END DO 
     
    230237            ! 
    231238            DO ji = 1, npti 
    232                ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    233                   &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    234                ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    235                   &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     239               ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     240                  &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     241               ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     242                  &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236243            END DO 
    237244            DO jk = 1, nlay_i-1 
    238245               DO ji = 1, npti 
    239                   ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    240                      &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
    241                      &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     246                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
     247                     &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
     248                     &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
    242249               END DO 
    243250            END DO 
    244251            ! 
    245252         ENDIF 
    246          ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )         
     253 
     254         ! Variable used after iterations 
     255         ! Value must be frozen after convergence for MPP independance reason 
     256         DO ji = 1, npti 
     257            IF ( .NOT. l_T_converged(ji) ) & 
     258               ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )         
     259         END DO 
    247260         ! 
    248261         !--- G(he) : enhancement of thermal conductivity in mono-category case 
     
    270283         !----------------- 
    271284         !--- Snow 
     285         ! Variable used after iterations 
     286         ! Value must be frozen after convergence for MPP independance reason 
    272287         DO jk = 0, nlay_s-1 
    273288            DO ji = 1, npti 
    274                zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     289               IF ( .NOT. l_T_converged(ji) ) & 
     290                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
    275291            END DO 
    276292         END DO 
    277293         DO ji = 1, npti   ! Snow-ice interface 
    278             zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
    279             IF( zfac > epsi10 ) THEN 
    280                zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
    281             ELSE 
    282                zkappa_s(ji,nlay_s) = 0._wp 
     294            IF ( .NOT. l_T_converged(ji) ) THEN 
     295               zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
     296               IF( zfac > epsi10 ) THEN 
     297                  zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
     298               ELSE 
     299                  zkappa_s(ji,nlay_s) = 0._wp 
     300               ENDIF 
    283301            ENDIF 
    284302         END DO 
    285303 
    286304         !--- Ice 
     305         ! Variable used after iterations 
     306         ! Value must be frozen after convergence for MPP independance reason 
    287307         DO jk = 0, nlay_i 
    288308            DO ji = 1, npti 
    289                zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
     309               IF ( .NOT. l_T_converged(ji) ) & 
     310                  zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
    290311            END DO 
    291312         END DO 
    292313         DO ji = 1, npti   ! Snow-ice interface 
    293             zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     314            IF ( .NOT. l_T_converged(ji) ) & 
     315               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    294316         END DO 
    295317         ! 
     
    326348            ! update of the non solar flux according to the update in T_su 
    327349            DO ji = 1, npti 
    328                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
     350               ! Variable used after iterations 
     351               ! Value must be frozen after convergence for MPP independance reason 
     352               IF ( .NOT. l_T_converged(ji) ) & 
     353                  qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    329354            END DO 
    330355 
     
    496521            ! ice temperatures 
    497522            DO ji = 1, npti 
    498                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     523               ! Variable used after iterations 
     524               ! Value must be frozen after convergence for MPP independance reason 
     525               IF ( .NOT. l_T_converged(ji) ) & 
     526                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    499527            END DO 
    500528 
     
    502530               DO ji = 1, npti 
    503531                  jk = jm - nlay_s - 1 
    504                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     532                  IF ( .NOT. l_T_converged(ji) ) & 
     533                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    505534               END DO 
    506535            END DO 
    507536 
    508537            DO ji = 1, npti 
    509                ! snow temperatures       
    510                IF( h_s_1d(ji) > 0._wp ) THEN 
    511                   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) 
    512                ENDIF 
    513                ! surface temperature 
    514                ztsub(ji) = t_su_1d(ji) 
    515                IF( t_su_1d(ji) < rt0 ) THEN 
    516                   t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    517                      &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     538               ! Variables used after iterations 
     539               ! Value must be frozen after convergence for MPP independance reason 
     540               IF ( .NOT. l_T_converged(ji) ) THEN 
     541                  ! snow temperatures       
     542                  IF( h_s_1d(ji) > 0._wp ) THEN 
     543                     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) 
     544                  ENDIF 
     545                  ! surface temperature 
     546                  ztsub(ji) = t_su_1d(ji) 
     547                  IF( t_su_1d(ji) < rt0 ) THEN 
     548                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     549                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     550                  ENDIF 
    518551               ENDIF 
    519552            END DO 
     
    524557            ! check that nowhere it has started to melt 
    525558            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    526             zdti_max = 0._wp 
    527             DO ji = 1, npti 
    528                t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    529                zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    530             END DO 
    531  
    532             DO jk = 1, nlay_s 
    533                DO ji = 1, npti 
    534                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    535                   zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    536                END DO 
    537             END DO 
    538  
    539             DO jk = 1, nlay_i 
    540                DO ji = 1, npti 
    541                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    542                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    543                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    544                END DO 
    545             END DO 
    546  
    547             ! Compute spatial maximum over all errors 
    548             ! note that this could be optimized substantially by iterating only the non-converging points 
    549             IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    550          ! 
     559 
     560            DO ji = 1, npti 
     561 
     562               zdti_max = 0._wp 
     563 
     564               IF ( .NOT. l_T_converged(ji) ) THEN 
     565                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
     566                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 
     567 
     568                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
     569                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
     570 
     571                  DO jk = 1, nlay_i 
     572                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
     573                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
     574                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     575                  END DO 
     576 
     577                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 
     578 
     579               ENDIF 
     580 
     581            END DO 
     582 
    551583         !----------------------------------------! 
    552584         !                                        ! 
     
    670702             
    671703            ! ice temperatures 
    672            DO ji = 1, npti 
    673                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     704            DO ji = 1, npti 
     705               ! Variable used after iterations 
     706               ! Value must be frozen after convergence for MPP independance reason 
     707               IF ( .NOT. l_T_converged(ji) ) & 
     708                  t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    674709            END DO 
    675710 
    676711            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    677712               DO ji = 1, npti 
    678                   jk = jm - nlay_s - 1 
    679                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     713                  IF ( .NOT. l_T_converged(ji) ) THEN 
     714                     jk = jm - nlay_s - 1 
     715                     t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     716                  ENDIF 
    680717               END DO 
    681718            END DO 
     
    683720            ! snow temperatures       
    684721            DO ji = 1, npti 
    685                IF( h_s_1d(ji) > 0._wp ) THEN 
    686                   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) 
     722               ! Variable used after iterations 
     723               ! Value must be frozen after convergence for MPP independance reason 
     724               IF ( .NOT. l_T_converged(ji) ) THEN 
     725                  IF( h_s_1d(ji) > 0._wp ) THEN 
     726                     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) 
     727                  ENDIF 
    687728               ENDIF 
    688729            END DO 
     
    693734            ! check that nowhere it has started to melt 
    694735            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    695             zdti_max = 0._wp 
    696  
    697             DO jk = 1, nlay_s 
    698                DO ji = 1, npti 
    699                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    700                   zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    701                END DO 
    702             END DO 
    703              
    704             DO jk = 1, nlay_i 
    705                DO ji = 1, npti 
    706                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    707                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    708                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    709                END DO 
    710             END DO 
    711  
    712             ! Compute spatial maximum over all errors 
    713             ! note that this could be optimized substantially by iterating only the non-converging points 
    714             IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     736 
     737            DO ji = 1, npti 
     738 
     739               zdti_max = 0._wp 
     740 
     741               IF ( .NOT. l_T_converged(ji) ) THEN 
     742                  ! t_s 
     743                  t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
     744                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
     745                  ! t_i 
     746                  DO jk = 0, nlay_i 
     747                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
     748                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
     749                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     750                  END DO 
     751 
     752                  IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 
     753 
     754               ENDIF 
     755 
     756            END DO 
    715757 
    716758         ENDIF ! k_jules 
Note: See TracChangeset for help on using the changeset viewer.