source: LMDZ4/branches/LMDZ4_V2_patch/libf/phylmd/calcul_simulISCCP.h @ 739

c
c $Header$
c
cIM 090704 BEG
c     nbapp_isccp=48
c     nbapp_isccp=8
c     nbapp_isccp=6
c     nbapp_isccp=4 !CPU < 30min pour  9pdt/jour
      nbapp_isccp=3 !CPU ??           10pdt/jour
c     nbapp_isccp=2
c     nbapp_isccp=1
      isccppas=NINT(86400./dtime/nbapp_isccp)
cIM 010904 BEG
cIM   IF (MOD(itap,isccppas).EQ.0) THEN
c      PRINT*,'itap,isccppas,xjour',itap,isccppas,xjour
c
cIM initialisation nbsunlit pour calculs simulateur ISCCP pdt la journee
c
      DO i=1, klon
       sunlit(i)=1 
       IF(rmu0(i).EQ.0.) sunlit(i)=0
       nbsunlit(1,i)=FLOAT(sunlit(i))
      ENDDO
c
cIM calcul tau, emissivite nuages convectifs
c
      convfra(:,:)=rnebcon(:,:)
      convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
cIM Amip2 beg
c
      CALL newmicro (paprs, pplay,ok_newmicro,
     .            t_seri, convliq, convfra, dtau_c, dem_c,
     .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
     .            flwp_c, fiwp_c, flwc_c, fiwc_c,
     e            ok_aie,
     e            sulfate, sulfate_pi,
     e            bl95_b0, bl95_b1,
     s            cldtaupi, re, fl)
c
cIM Amip2 end

c
cIM calcul tau, emissivite nuages startiformes
c
cIM Amip2 beg
c
      CALL newmicro (paprs, pplay,ok_newmicro,
     .            t_seri, cldliq, cldfra, dtau_s, dem_s,
     .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
     .            flwp_s, fiwp_s, flwc_s, fiwc_s,
     e            ok_aie,
     e            sulfate, sulfate_pi,
     e            bl95_b0, bl95_b1,
     s            cldtaupi, re, fl)
c
cIM Amip2 end
c
      cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
c
cIM inversion des niveaux de pression ==> de haut en bas
c
      CALL haut2bas(klon, klev, pplay, pfull)
      CALL haut2bas(klon, klev, q_seri, qv)
      CALL haut2bas(klon, klev, cldtot, cc)
      CALL haut2bas(klon, klev, rnebcon, conv)
      CALL haut2bas(klon, klev, dtau_s, dtau_sH2B)
      CALL haut2bas(klon, klev, dtau_c, dtau_cH2B)
      CALL haut2bas(klon, klev, t_seri, at)
      CALL haut2bas(klon, klev, dem_s, dem_sH2B)
      CALL haut2bas(klon, klev, dem_c, dem_cH2B)
      CALL haut2bas(klon, klevp1, paprs, phalf)
c
cIM lecture invtau, tautab des fichiers formattes
c
      IF (debut) THEN
c     open(99,file='tautab.bin',access='sequential',
c    $     form='unformatted',status='old')
c     read(99) tautab
c
      open(99,file='tautab.formatted', FORM='FORMATTED')
      read(99,'(f30.20)') tautab
      close(99)
c
      open(99,file='invtau.formatted',form='FORMATTED')
      read(99,'(i10)') invtau
      close(99)
c
cIM: calcul coordonnees regions pour statistiques distribution 
cIM: nuages en ftion du regime dynamique pour regions oceaniques
c
      IF (ok_regdyn) THEN !histREGDYN
c
#include "ini_coord_REGDYN.h"
c
      ENDIF !ok_regdyn
c 
      ENDIF !debut
c
cIM: initialisation de seed 
c       DO i=1, klon
c         seed(i)=i+100
c       ENDDO
c
        DO i=1, klon
c
         aa=ABS(paprs(i,2)-NINT(paprs(i,2)))
         seed_re(i)=1000.*aa+1.
         seed(i)=NINT(seed_re(i))
c
         IF(seed(i).LT.50) THEN
c          print*,'seed<50 avant i seed itap paprs',i,
c    .     seed(i),itap,paprs(i,2)
           seed(i)=50+seed(i)+i+itap
           seed_old(i)=seed(i)
c
           IF(itap.GT.1) then
            IF(seed(i).EQ.seed_old(i)) THEN
             seed(i)=seed(i)+10
             seed_old(i)=seed(i)
            ENDIF   
           ENDIF
c
c          print*,'seed<50 apres i seed itap paprs',i,
c    .     seed(i),itap,paprs(i,2)
c
          ELSE IF(seed(i).EQ.0) THEN
           print*,'seed=0 i paprs aa seed_re',
     .     i,paprs(i,2),aa,seed_re(i)
           STOP    
          ELSE IF(seed(i).LT.0) THEN
           print*,'seed < 0, i seed itap paprs',i,
     .     seed(i),itap,paprs(i,2)
           STOP    
          ENDIF
c
        ENDDO
c     
cIM: pas de debug, debugcol
      debug=0
      debugcol=0
c
cIM o500 ==> distribution nuage ftion du regime dynamique (vit. verticale a 500 hPa)
c
        DO k=1, klevm1
        kp1=k+1
c       PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1)
        if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN
         DO i=1, klon
          o500(i)=omega(i,k)*RDAY/100.
c         if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1)
         ENDDO
         GOTO 1000
        endif 
1000  continue
      ENDDO
c
cIM recalcule les nuages vus par satellite, via le simulateur ISCCP
c
      CALL ISCCP_CLOUD_TYPES(
     &     debug,
     &     debugcol,
cIM 300704    &     itap, debut,
cIM 300604 klon !BAD
     &     klon,
     &     sunlit,
     &     klev,
     &     ncol,
     &     seed,
     &     pfull,
     &     phalf,
     &     qv, cc, conv, dtau_sH2B, dtau_cH2B,
     &     top_height,
     &     overlap,
     &     tautab,
     &     invtau,
     &     ztsol,
     &     emsfc_lw,
     &     at, dem_sH2B, dem_cH2B,
     &     fq_isccp,
     &     totalcldarea,
     &     meanptop,
     &     meantaucld,
     &     boxtau,
     &     boxptop)
c
c calcul regime dynamique sur les regions fixees
c
       IF (ok_regdyn) THEN !histREGDYN
c
#include "calcul_REGDYN.h"
c
       ENDIF !(ok_regdyn) THEN !histREGDYN
cIM    ENDIF !(MOD(itaprad,radpas).EQ.0) THEN
cIM 010904 END