Ignore:
Timestamp:
Feb 3, 2009, 11:26:58 AM (16 years ago)
Author:
yann meurdesoif
Message:

Optimisation Othman Bouizi : newmicro

YM

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F

    r766 r1092  
    110110      REAL zclear(klon)
    111111      REAL zcloud(klon)
     112
     113c **************************
     114c *                        *
     115c * DEBUT PARTIE OPTIMISEE *
     116c *                        *
     117c **************************
     118
     119      REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
     120      REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
     121
    112122c
    113123c Calculer l'epaisseur optique et l'emmissivite des nuages
    114124c
    115 cIM inversion des DO
    116       DO i = 1, klon
    117        xflwp(i)=0.
    118        xfiwp(i)=0.
    119 cccccccccccc!CDIR NOVECTOR
     125c     IM inversion des DO
     126      xflwp = 0.d0
     127      xfiwp = 0.d0
     128      xflwc = 0.d0
     129      xfiwc = 0.d0
     130
    120131      DO k = 1, klev
    121 c
    122        xflwc(i,k)=0.
    123        xfiwc(i,k)=0.
    124 c
    125          rad_chaud = rad_chau1
    126          IF (k.LE.3) rad_chaud = rad_chau2
    127          pclc(i,k) = MAX(pclc(i,k), seuil_neb)
    128          zflwp(i) = 1000.*pqlwp(i,k)/RG/pclc(i,k)
    129      .          *(paprs(i,k)-paprs(i,k+1))
    130          zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
    131          zfice = MIN(MAX(zfice,0.0),1.0)
    132          zfice = zfice**nexpo
    133          radius = rad_chaud * (1.-zfice) + rad_froid * zfice
    134          coef = coef_chau * (1.-zfice) + coef_froi * zfice
    135          pcltau(i,k) = 3.0/2.0 * zflwp(i) / radius
    136          pclemi(i,k) = 1.0 - EXP( - coef * zflwp(i))
    137 
    138          if (ok_newmicro) then
    139 
    140 c -- liquid/ice cloud water paths:
    141 
    142          zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
    143          zfice = MIN(MAX(zfice,0.0),1.0)
    144 
    145          zflwp(i) = 1000.*(1.-zfice)*pqlwp(i,k)/pclc(i,k)
    146      :          *(paprs(i,k)-paprs(i,k+1))/RG
    147          zfiwp(i) = 1000.*zfice*pqlwp(i,k)/pclc(i,k)
    148      :          *(paprs(i,k)-paprs(i,k+1))/RG
    149 
    150          xflwp(i) = xflwp(i)+ (1.-zfice)*pqlwp(i,k)
    151      :          *(paprs(i,k)-paprs(i,k+1))/RG
    152          xfiwp(i) = xfiwp(i)+ zfice*pqlwp(i,k)
    153      :          *(paprs(i,k)-paprs(i,k+1))/RG
    154 
    155 cIM Total Liquid/Ice water content
    156          xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)
    157          xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)
    158 cIM In-Cloud Liquid/Ice water content
    159 c        xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
    160 c        xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
    161 
    162 c -- effective cloud droplet radius (microns):
    163 
    164 c for liquid water clouds:
     132         DO i = 1, klon
     133            diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
     134         ENDDO
     135      ENDDO
     136
     137      IF (ok_newmicro) THEN
     138
     139
     140         DO k = 1, klev
     141            DO i = 1, klon
     142               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
     143               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
     144c     IM Total Liquid/Ice water content                                   
     145               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
     146               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
     147c     IM In-Cloud Liquid/Ice water content
     148c     xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
     149c     xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
     150            ENDDO
     151         ENDDO
     152
    165153         IF (ok_aie) THEN
    166             ! Formula "D" of Boucher and Lohmann, Tellus, 1995
    167             !             
    168             cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
    169      .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    170             ! Cloud droplet number concentration (CDNC) is restricted
    171             ! to be within [20, 1000 cm^3]
    172             !
    173             cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
    174             !
    175             !
    176             cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
    177      .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    178             cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
    179             !           
    180             !
    181             ! air density: pplay(i,k) / (RD * zT(i,k))
    182             ! factor 1.1: derive effective radius from volume-mean radius
    183             ! factor 1000 is the water density
    184             ! _chaud means that this is the CDR for liquid water clouds
    185             !
    186             rad_chaud =
    187      .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
    188      .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
    189             !
    190             ! Convert to um. CDR shall be at least 3 um.
    191             !
    192 c           rad_chaud = MAX(rad_chaud*1.e6, 3.)
    193             rad_chaud = MAX(rad_chaud*1.e6, 5.)
    194            
    195             ! Pre-industrial cloud opt thickness
    196             !
    197             ! "radius" is calculated as rad_chaud above (plus the
    198             ! ice cloud contribution) but using cdnc_pi instead of
    199             ! cdnc.
    200             radius =
    201      .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
    202      .               / (4./3. * RPI * 1000. * cdnc_pi(i,k)) )**(1./3.)
    203             radius = MAX(radius*1.e6, 5.)
    204            
    205             tc = t(i,k)-273.15
    206             rei = 0.71*tc + 61.29
    207             if (tc.le.-81.4) rei = 3.5
    208             if (zflwp(i).eq.0.) radius = 1.
    209             if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
    210             cldtaupi(i,k) = 3.0/2.0 * zflwp(i) / radius
    211      .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
    212          ENDIF                  ! ok_aie
    213          ! For output diagnostics
    214          !
    215          ! Cloud droplet effective radius [um]
    216          !
    217          ! we multiply here with f * xl (fraction of liquid water
    218          ! clouds in the grid cell) to avoid problems in the
    219          ! averaging of the output.
    220          ! In the output of IOIPSL, derive the real cloud droplet
    221          ! effective radius as re/fl
    222          !
    223          fl(i,k) = pclc(i,k)*(1.-zfice)           
    224          re(i,k) = rad_chaud*fl(i,k)
    225            
    226 c-jq end         
     154            DO k = 1, klev
     155               DO i = 1, klon
     156                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
     157                                !             
     158                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
     159     &                 log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     160                                ! Cloud droplet number concentration (CDNC) is restricted
     161                                ! to be within [20, 1000 cm^3]
     162                                !
     163                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
     164                                !
     165                                !
     166                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
     167     &                 log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     168                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
     169               ENDDO
     170            ENDDO
     171            DO k = 1, klev
     172               DO i = 1, klon
     173!                  rad_chaud_tab(i,k) =
     174!     &                 MAX(1.1e6
     175!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
     176!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
     177                  rad_chaud_tab(i,k) =
     178     &                 1.1
     179     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
     180     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
     181                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
     182               ENDDO           
     183            ENDDO
     184         ELSE
     185            DO k = 1, MIN(3,klev)
     186               DO i = 1, klon
     187                  rad_chaud_tab(i,k) = rad_chau2
     188               ENDDO           
     189            ENDDO
     190            DO k = MIN(3,klev)+1, klev
     191               DO i = 1, klon
     192                  rad_chaud_tab(i,k) = rad_chau1
     193               ENDDO           
     194            ENDDO
     195
     196         ENDIF
    227197         
    228          rel = rad_chaud
    229 c for ice clouds: as a function of the ambiant temperature
    230 c [formula used by Iacobellis and Somerville (2000), with an
    231 c asymptotical value of 3.5 microns at T<-81.4 C added to be
    232 c consistent with observations of Heymsfield et al. 1986]:
    233          tc = t(i,k)-273.15
    234          rei = 0.71*tc + 61.29
    235          if (tc.le.-81.4) rei = 3.5
    236 
    237 c -- cloud optical thickness :
    238 
    239 c [for liquid clouds, traditional formula,
    240 c  for ice clouds, Ebert & Curry (1992)]
    241 
    242          if (zflwp(i).eq.0.) rel = 1.
    243          if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
    244          pcltau(i,k) = 3.0/2.0 * ( zflwp(i)/rel )
    245      .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
    246 
    247 c -- cloud infrared emissivity:
    248 
    249 c [the broadband infrared absorption coefficient is parameterized
    250 c  as a function of the effective cld droplet radius]
    251 
    252 c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
    253          k_ice = k_ice0 + 1.0/rei
    254 
    255          pclemi(i,k) = 1.0
    256      .      - EXP( - coef_chau*zflwp(i) - DF*k_ice*zfiwp(i) )
    257 
    258          endif ! ok_newmicro
    259 
    260          lo = (pclc(i,k) .LE. seuil_neb)
    261          IF (lo) pclc(i,k) = 0.0
    262          IF (lo) pcltau(i,k) = 0.0
    263          IF (lo) pclemi(i,k) = 0.0
    264          
    265          IF (lo) cldtaupi(i,k) = 0.0
    266          IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)           
    267       ENDDO
    268       ENDDO
    269 ccc      DO k = 1, klev
    270 ccc      DO i = 1, klon
    271 ccc         t(i,k) = t(i,k)
    272 ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
    273 ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
    274 ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
    275 ccc     .          /(rg*pclc(i,k))
    276 ccc         zradef = 10.0 + (1.-sigs(k))*45.0
    277 ccc         pcltau(i,k) = 1.5 * zflwp / zradef
    278 ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
    279 ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
    280 ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
    281 ccc         if (.NOT.lo) pclc(i,k) = 0.0
    282 ccc         if (.NOT.lo) pcltau(i,k) = 0.0
    283 ccc         if (.NOT.lo) pclemi(i,k) = 0.0
    284 ccc      ENDDO
    285 ccc      ENDDO
    286 cccccc      print*, 'pas de nuage dans le rayonnement'
    287 cccccc      DO k = 1, klev
    288 cccccc      DO i = 1, klon
    289 cccccc         pclc(i,k) = 0.0
    290 cccccc         pcltau(i,k) = 0.0
    291 cccccc         pclemi(i,k) = 0.0
    292 cccccc      ENDDO
    293 cccccc      ENDDO
    294 C
    295 C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
    296 C
    297 cIM cf. CR:test: calcul prenant ou non en compte le recouvrement
    298 cinitialisations
     198         DO k = 1, klev
     199!            IF(.not.ok_aie) THEN
     200            rad_chaud = rad_chau1
     201            IF (k.LE.3) rad_chaud = rad_chau2
     202!            ENDIF
     203            DO i = 1, klon
     204               IF (pclc(i,k) .LE. seuil_neb) THEN
     205               
     206c     -- effective cloud droplet radius (microns):
     207               
     208c     for liquid water clouds:
     209                                ! For output diagnostics
     210                                !
     211                                ! Cloud droplet effective radius [um]
     212                                !
     213                                ! we multiply here with f * xl (fraction of liquid water
     214                                ! clouds in the grid cell) to avoid problems in the
     215                                ! averaging of the output.
     216                                ! In the output of IOIPSL, derive the real cloud droplet
     217                                ! effective radius as re/fl
     218                                !
     219                                   
     220                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))           
     221                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
     222                 
     223                  pclc(i,k) = 0.0
     224                  pcltau(i,k) = 0.0
     225                  pclemi(i,k) = 0.0
     226                  cldtaupi(i,k) = 0.0                 
     227               ELSE
     228
     229c     -- liquid/ice cloud water paths:
     230                 
     231                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
     232     &                 *diff_paprs(i,k)
     233                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
     234     &                 *diff_paprs(i,k)
     235                 
     236c     -- effective cloud droplet radius (microns):
     237               
     238c     for liquid water clouds:
     239                                   
     240                  IF (ok_aie) THEN
     241                     radius =
     242     &                    1.1
     243     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
     244     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
     245                     radius = MAX(radius*1e6, 5.)
     246                 
     247                     tc = t(i,k)-273.15
     248                     rei = 0.71*tc + 61.29
     249                     if (tc.le.-81.4) rei = 3.5
     250                     if (zflwp_var.eq.0.) radius = 1.
     251                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
     252                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
     253     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
     254                  ENDIF         ! ok_aie
     255                                ! For output diagnostics
     256                                !
     257                                ! Cloud droplet effective radius [um]
     258                                !
     259                                ! we multiply here with f * xl (fraction of liquid water
     260                                ! clouds in the grid cell) to avoid problems in the
     261                                ! averaging of the output.
     262                                ! In the output of IOIPSL, derive the real cloud droplet
     263                                ! effective radius as re/fl
     264                                !
     265 
     266                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))           
     267                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
     268                 
     269                  rel = rad_chaud_tab(i,k)
     270c     for ice clouds: as a function of the ambiant temperature
     271c     [formula used by Iacobellis and Somerville (2000), with an
     272c     asymptotical value of 3.5 microns at T<-81.4 C added to be
     273c     consistent with observations of Heymsfield et al. 1986]:
     274                  tc = t(i,k)-273.15
     275                  rei = 0.71*tc + 61.29
     276                  if (tc.le.-81.4) rei = 3.5
     277c     -- cloud optical thickness :
     278               
     279c     [for liquid clouds, traditional formula,
     280c     for ice clouds, Ebert & Curry (1992)]
     281                 
     282                  if (zflwp_var.eq.0.) rel = 1.
     283                  if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
     284                  pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
     285     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
     286c     -- cloud infrared emissivity:
     287               
     288c     [the broadband infrared absorption coefficient is parameterized
     289c     as a function of the effective cld droplet radius]
     290               
     291c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
     292                  k_ice = k_ice0 + 1.0/rei
     293                 
     294                  pclemi(i,k) = 1.0
     295     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
     296
     297               ENDIF
     298               
     299            ENDDO
     300         ENDDO
     301
     302         DO k = 1, klev
     303            DO i = 1, klon
     304               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
     305               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
     306            ENDDO
     307         ENDDO
     308
     309      ELSE
     310         DO k = 1, klev
     311            rad_chaud = rad_chau1
     312            IF (k.LE.3) rad_chaud = rad_chau2
     313            DO i = 1, klon
     314                             
     315               IF (pclc(i,k) .LE. seuil_neb) THEN
     316
     317                  pclc(i,k) = 0.0
     318                  pcltau(i,k) = 0.0
     319                  pclemi(i,k) = 0.0
     320                  cldtaupi(i,k) = 0.0
     321
     322               ELSE
     323
     324                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
     325     &                 /pclc(i,k)
     326                 
     327                  zfice1 = MIN(
     328     &                 MAX( 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
     329     &                 ,0.0),1.0)**nexpo
     330                 
     331                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
     332                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
     333
     334                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
     335                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
     336
     337               ENDIF
     338                             
     339            ENDDO
     340         ENDDO
     341      ENDIF
     342     
     343      IF (.NOT.ok_aie) THEN
     344         DO k = 1, klev
     345            DO i = 1, klon
     346               cldtaupi(i,k)=pcltau(i,k)
     347            ENDDO
     348         ENDDO               
     349      ENDIF
     350
     351ccc   DO k = 1, klev
     352ccc   DO i = 1, klon
     353ccc   t(i,k) = t(i,k)
     354ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
     355ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
     356ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
     357ccc   .          /(rg*pclc(i,k))
     358ccc   zradef = 10.0 + (1.-sigs(k))*45.0
     359ccc   pcltau(i,k) = 1.5 * zflwp / zradef
     360ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
     361ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
     362ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
     363ccc   if (.NOT.lo) pclc(i,k) = 0.0
     364ccc   if (.NOT.lo) pcltau(i,k) = 0.0
     365ccc   if (.NOT.lo) pclemi(i,k) = 0.0
     366ccc   ENDDO
     367ccc   ENDDO
     368ccccc print*, 'pas de nuage dans le rayonnement'
     369ccccc DO k = 1, klev
     370ccccc DO i = 1, klon
     371ccccc pclc(i,k) = 0.0
     372ccccc pcltau(i,k) = 0.0
     373ccccc pclemi(i,k) = 0.0
     374ccccc ENDDO
     375ccccc ENDDO
     376C     
     377C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
     378C     
     379c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
     380c     initialisations
    299381      DO i=1,klon
    300382         zclear(i)=1.
     
    308390cIM cf CR DO k=1,klev
    309391      DO k = klev, 1, -1
    310       DO i = 1, klon
    311          pctlwp(i) = pctlwp(i)
    312      .             + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
    313 cIM cf. CR
    314             IF (NOVLP.EQ.1) THEN
     392         DO i = 1, klon
     393            pctlwp(i) = pctlwp(i)
     394     &           + pqlwp(i,k)*diff_paprs(i,k)
     395         ENDDO
     396      ENDDO
     397c     IM cf. CR
     398      IF (NOVLP.EQ.1) THEN
     399         DO k = klev, 1, -1
     400            DO i = 1, klon
    315401               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
    316      s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
     402     &              /(1.-MIN(zcloud(i),1.-ZEPSEC))
    317403               pct(i)=1.-zclear(i)
    318                if (pplay(i,k).LE.cetahb*paprs(i,1)) then
     404               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
    319405                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
    320      s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
    321                else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
    322      .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
     406     &                 /(1.-MIN(zcloud(i),1.-ZEPSEC))
     407               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
     408     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
    323409                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
    324      s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
    325                else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
     410     &                 /(1.-MIN(zcloud(i),1.-ZEPSEC))
     411               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
    326412                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
    327      s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
     413     &                 /(1.-MIN(zcloud(i),1.-ZEPSEC))
    328414               endif
    329415               zcloud(i)=pclc(i,k)
    330             ELSE IF (NOVLP.EQ.2) THEN
     416            ENDDO
     417         ENDDO
     418      ELSE IF (NOVLP.EQ.2) THEN
     419         DO k = klev, 1, -1
     420            DO i = 1, klon
    331421               zcloud(i)=MAX(pclc(i,k),zcloud(i))
    332422               pct(i)=zcloud(i)
    333                if (pplay(i,k).LE.cetahb*paprs(i,1)) then
     423               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
    334424                  pch(i) = MIN(pclc(i,k),pch(i))
    335                else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
    336      .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
     425               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
     426     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
    337427                  pcm(i) = MIN(pclc(i,k),pcm(i))
    338                else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
     428               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
    339429                  pcl(i) = MIN(pclc(i,k),pcl(i))
    340430               endif
    341             ELSE IF (NOVLP.EQ.3) THEN
     431            ENDDO
     432         ENDDO
     433      ELSE IF (NOVLP.EQ.3) THEN
     434         DO k = klev, 1, -1
     435            DO i = 1, klon
    342436               zclear(i)=zclear(i)*(1.-pclc(i,k))
    343437               pct(i)=1-zclear(i)
    344                if (pplay(i,k).LE.cetahb*paprs(i,1)) then
    345                pch(i) = pch(i)*(1.0-pclc(i,k))
    346                else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
    347      .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
    348                pcm(i) = pcm(i)*(1.0-pclc(i,k))
    349                else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
    350                pcl(i) = pcl(i)*(1.0-pclc(i,k))
     438               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
     439                  pch(i) = pch(i)*(1.0-pclc(i,k))
     440               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
     441     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
     442                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
     443               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
     444                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
    351445               endif
    352             ENDIF
    353          ENDDO
    354       ENDDO
    355 C
     446            ENDDO
     447         ENDDO
     448      ENDIF
     449     
     450C     
    356451      DO i = 1, klon
    357 cIM cf. CR          pct(i)=1.-pct(i)
     452c     IM cf. CR          pct(i)=1.-pct(i)
    358453         pch(i)=1.-pch(i)
    359454         pcm(i)=1.-pcm(i)
    360455         pcl(i)=1.-pcl(i)
    361456      ENDDO
     457     
    362458C
    363459      RETURN
Note: See TracChangeset for help on using the changeset viewer.