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 9454 – NEMO

Changeset 9454


Ignore:
Timestamp:
2018-04-03T17:33:37+02:00 (6 years ago)
Author:
clem
Message:

make sure that agrif is working with sea ice in any conditions

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90

    r9167 r9454  
    2424   USE ice 
    2525   USE agrif_ice 
     26   USE phycst , ONLY: rt0 
    2627    
    2728   IMPLICIT NONE 
     
    183184      ELSE               ! child grid 
    184185         ! 
    185          IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used 
     186!         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used 
    186187            ! 
    187188            jm = 1 
     
    214215            END DO 
    215216            ! 
    216          ELSE                          ! ==> complex interpolation (only one ghost cell available) 
    217             !! Use a more complex interpolation since we mix solutions over a couple of grid points 
    218             !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...) 
    219             !! clem: for some reason (I don't know why), the following lines do not work  
    220             !!       with mpp (or in realistic configurations?). It makes the model crash 
    221             !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly) 
    222             ! record ztab 
    223             jm = 1 
    224             DO jl = 1, jpl 
    225                ztab(:,:,jm  ) = a_i (:,:,jl) 
    226                ztab(:,:,jm+1) = v_i (:,:,jl) 
    227                ztab(:,:,jm+2) = v_s (:,:,jl) 
    228                ztab(:,:,jm+3) = sv_i(:,:,jl) 
    229                ztab(:,:,jm+4) = oa_i(:,:,jl) 
    230                ztab(:,:,jm+5) = a_ip(:,:,jl) 
    231                ztab(:,:,jm+6) = v_ip(:,:,jl) 
    232                ztab(:,:,jm+7) = t_su(:,:,jl) 
    233                jm = jm + 8 
    234                DO jk = 1, nlay_s 
    235                   ztab(:,:,jm) = e_s(:,:,jk,jl) 
    236                   jm = jm + 1 
    237                END DO 
    238                DO jk = 1, nlay_i 
    239                   ztab(:,:,jm) = e_i(:,:,jk,jl) 
    240                   jm = jm + 1 
    241                END DO 
    242                ! 
    243             END DO 
    244             ! 
    245             ! borders of the domain 
    246             western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2) 
    247             southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2) 
    248             ! 
    249             ! spatial smoothing 
    250             zrhox = Agrif_Rhox() 
    251             z1 =      ( zrhox - 1. ) * 0.5  
    252             z3 =      ( zrhox - 1. ) / ( zrhox + 1. ) 
    253             z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
    254             z7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
    255             z2 = 1. - z1 
    256             z4 = 1. - z3 
    257             z5 = 1. - z6 - z7 
    258             ! 
    259             ! Remove corners 
    260             imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2 
    261             IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3 
    262             IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2 
    263             IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3 
    264             IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2 
    265  
    266             ! smoothed fields 
    267             IF( eastern_side ) THEN 
    268                ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:) 
    269                DO jj = jmin, jmax 
    270                   rswitch = 0. 
    271                   IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1. 
    272                   ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  & 
    273                      &                +      umask(nlci-2,jj,1)   *  & 
    274                      &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  & 
    275                      &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) ) 
    276                   ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1) 
    277                END DO 
    278             ENDIF 
    279             !  
    280             IF( northern_side ) THEN 
    281                ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:) 
    282                DO ji = imin, imax 
    283                   rswitch = 0. 
    284                   IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1. 
    285                   ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  & 
    286                      &                +      vmask(ji,nlcj-2,1)   *  & 
    287                      &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) & 
    288                      &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) ) 
    289                   ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1) 
    290                END DO 
    291             END IF 
    292             ! 
    293             IF( western_side) THEN 
    294                ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:) 
    295                DO jj = jmin, jmax 
    296                   rswitch = 0. 
    297                   IF( u_ice(2,jj) < 0._wp ) rswitch = 1. 
    298                   ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  & 
    299                      &           +      umask(2,jj,1)   *   & 
    300                      &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) & 
    301                      &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) ) 
    302                   ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1) 
    303                END DO 
    304             ENDIF 
    305             ! 
    306             IF( southern_side ) THEN 
    307                ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:) 
    308                DO ji = imin, imax 
    309                   rswitch = 0. 
    310                   IF( v_ice(ji,2) < 0._wp ) rswitch = 1. 
    311                   ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  & 
    312                      &           +      vmask(ji,2,1)   *  & 
    313                      &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) & 
    314                      &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) ) 
    315                   ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1) 
    316                END DO 
    317             END IF 
    318             ! 
    319             ! Treatment of corners 
    320             IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south 
    321             IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north 
    322             IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south 
    323             IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north 
    324              
    325             ! retrieve ice tracers 
    326             jm = 1 
    327             DO jl = 1, jpl 
    328                ! 
    329                DO jj = j1, j2 
    330                   DO ji = i1, i2 
    331                      a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1) 
    332                      v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1) 
    333                      v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1) 
    334                      sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1) 
    335                      oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1) 
    336                      a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1) 
    337                      v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1) 
    338                      t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1) 
    339                   END DO 
    340                END DO 
    341                jm = jm + 8 
    342                ! 
    343                DO jk = 1, nlay_s 
    344                   e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
    345                   jm = jm + 1 
    346                END DO 
    347                ! 
    348                DO jk = 1, nlay_i 
    349                   e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
    350                   jm = jm + 1 
    351                END DO 
    352                ! 
    353             END DO 
    354            
    355          ENDIF  ! nbghostcells=1 
     217!!==> clem: this interpolation does not work because it creates negative values, due 
     218!           to negative coefficients when mixing points (for ex. z7) 
     219!! 
     220!         ELSE                          ! ==> complex interpolation (only one ghost cell available) 
     221!            !! Use a more complex interpolation since we mix solutions over a couple of grid points 
     222!            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...) 
     223!            !! clem: for some reason (I don't know why), the following lines do not work  
     224!            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly) 
     225!            ! record ztab 
     226!            jm = 1 
     227!            DO jl = 1, jpl 
     228!               ztab(:,:,jm  ) = a_i (:,:,jl) 
     229!               ztab(:,:,jm+1) = v_i (:,:,jl) 
     230!               ztab(:,:,jm+2) = v_s (:,:,jl) 
     231!               ztab(:,:,jm+3) = sv_i(:,:,jl) 
     232!               ztab(:,:,jm+4) = oa_i(:,:,jl) 
     233!               ztab(:,:,jm+5) = a_ip(:,:,jl) 
     234!               ztab(:,:,jm+6) = v_ip(:,:,jl) 
     235!               ztab(:,:,jm+7) = t_su(:,:,jl) 
     236!               jm = jm + 8 
     237!               DO jk = 1, nlay_s 
     238!                  ztab(:,:,jm) = e_s(:,:,jk,jl) 
     239!                  jm = jm + 1 
     240!               END DO 
     241!               DO jk = 1, nlay_i 
     242!                  ztab(:,:,jm) = e_i(:,:,jk,jl) 
     243!                  jm = jm + 1 
     244!               END DO 
     245!               ! 
     246!            END DO 
     247!            ! 
     248!            ! borders of the domain 
     249!            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2) 
     250!            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2) 
     251!            ! 
     252!            ! spatial smoothing 
     253!            zrhox = Agrif_Rhox() 
     254!            z1 =      ( zrhox - 1. ) * 0.5  
     255!            z3 =      ( zrhox - 1. ) / ( zrhox + 1. ) 
     256!            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
     257!            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
     258!            z2 = 1. - z1 
     259!            z4 = 1. - z3 
     260!            z5 = 1. - z6 - z7 
     261!            ! 
     262!            ! Remove corners 
     263!            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2 
     264!            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3 
     265!            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2 
     266!            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3 
     267!            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2 
     268! 
     269!            ! smoothed fields 
     270!            IF( eastern_side ) THEN 
     271!               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:) 
     272!               DO jj = jmin, jmax 
     273!                  rswitch = 0. 
     274!                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1. 
     275!                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  & 
     276!                     &                +      umask(nlci-2,jj,1)   *  & 
     277!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  & 
     278!                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) ) 
     279!                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1) 
     280!               END DO 
     281!            ENDIF 
     282!            !  
     283!            IF( northern_side ) THEN 
     284!               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:) 
     285!               DO ji = imin, imax 
     286!                  rswitch = 0. 
     287!                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1. 
     288!                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  & 
     289!                     &                +      vmask(ji,nlcj-2,1)   *  & 
     290!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) & 
     291!                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) ) 
     292!                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1) 
     293!               END DO 
     294!            END IF 
     295!            ! 
     296!            IF( western_side) THEN 
     297!               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:) 
     298!               DO jj = jmin, jmax 
     299!                  rswitch = 0. 
     300!                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1. 
     301!                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  & 
     302!                     &           +      umask(2,jj,1)   *   & 
     303!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) & 
     304!                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) ) 
     305!                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1) 
     306!               END DO 
     307!            ENDIF 
     308!            ! 
     309!            IF( southern_side ) THEN 
     310!               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:) 
     311!               DO ji = imin, imax 
     312!                  rswitch = 0. 
     313!                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1. 
     314!                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  & 
     315!                     &           +      vmask(ji,2,1)   *  & 
     316!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) & 
     317!                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) ) 
     318!                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1) 
     319!               END DO 
     320!            END IF 
     321!            ! 
     322!            ! Treatment of corners 
     323!            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south 
     324!            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north 
     325!            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south 
     326!            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north 
     327!             
     328!            ! retrieve ice tracers 
     329!            jm = 1 
     330!            DO jl = 1, jpl 
     331!               ! 
     332!               DO jj = j1, j2 
     333!                  DO ji = i1, i2 
     334!                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1) 
     335!                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1) 
     336!                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1) 
     337!                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1) 
     338!                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1) 
     339!                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1) 
     340!                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1) 
     341!                     t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1) 
     342!                  END DO 
     343!               END DO 
     344!               jm = jm + 8 
     345!               ! 
     346!               DO jk = 1, nlay_s 
     347!                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
     348!                  jm = jm + 1 
     349!               END DO 
     350!               ! 
     351!               DO jk = 1, nlay_i 
     352!                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) 
     353!                  jm = jm + 1 
     354!               END DO 
     355!               ! 
     356!            END DO 
     357!           
     358!         ENDIF  ! nbghostcells=1 
    356359          
    357360         ! integrated values 
     
    361364         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 ) 
    362365         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 ) 
     366 
     367         at_ip(i1:i2,j1:j2) = SUM( a_ip(i1:i2,j1:j2,:), dim=3 ) ! melt ponds 
     368         vt_ip(i1:i2,j1:j2) = SUM( v_ip(i1:i2,j1:j2,:), dim=3 ) 
     369         ! 
     370         ato_i(i1:i2,j1:j2) = 1._wp - at_i(i1:i2,j1:j2)         ! open water fraction   
     371 
     372         DO jl = 1, jpl 
     373            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   t_su(i1:i2,j1:j2,jl) = rt0   ! to avoid a division by 0 in sbcblk.F90 
     374         END DO 
    363375         ! 
    364376      ENDIF 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90

    r9167 r9454  
    2727   USE ice 
    2828   USE agrif_ice  
     29   USE phycst , ONLY: rt0 
    2930 
    3031   IMPLICIT NONE 
     
    158159         ! 
    159160         ato_i(i1:i2,j1:j2) = 1._wp - at_i(i1:i2,j1:j2)         ! open water fraction   
     161 
     162         DO jl = 1, jpl 
     163            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   t_su(i1:i2,j1:j2,jl) = rt0   ! to avoid a division by 0 in sbcblk.F90 
     164         END DO 
    160165          
    161166      ENDIF 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r9098 r9454  
    1010   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines 
    1111   !!---------------------------------------------------------------------- 
    12  
     12#if defined key_agrif 
     13!DIR$ OPTIMIZE:1 
     14#endif 
    1315   !!---------------------------------------------------------------------- 
    1416   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r9190 r9454  
    825825      !! 
    826826      INTEGER  ::   ji, jj, jl               ! dummy loop indices 
    827       REAL(wp) ::   zst2, zst3               ! local variable 
     827      REAL(wp) ::   zst3                     ! local variable 
    828828      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    829829      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
    830830      REAL(wp) ::   zfr1, zfr2               ! local variables 
     831      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
    831832      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
    832833      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice 
     
    843844      ! 
    844845      zztmp = 1. / ( 1. - albo ) 
     846      WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
     847      ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp 
     848      END WHERE 
    845849      !                                     ! ========================== ! 
    846850      DO jl = 1, jpl                        !  Loop over ice categories  ! 
     
    851855               !      I   Radiative FLUXES   ! 
    852856               ! ----------------------------! 
    853                zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
    854                zst3 = ptsu(ji,jj,jl) * zst2 
     857               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
    855858               ! Short Wave (sw) 
    856859               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
     
    869872               ! Latent Heat 
    870873               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    871                   &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     874                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    872875               ! Latent heat sensitivity for ice (Dqla/Dt) 
    873876               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    874                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl)) 
     877                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     878                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl)) 
    875879               ELSE 
    876880                  dqla_ice(ji,jj,jl) = 0._wp 
Note: See TracChangeset for help on using the changeset viewer.