/[lmdze]/trunk/phylmd/physiq.f
ViewVC logotype

Diff of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 68 by guez, Wed Nov 14 16:59:30 2012 UTC revision 71 by guez, Mon Jul 8 18:12:18 2013 UTC
# Line 56  contains Line 56  contains
56      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
57      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
58      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
59        use readsulfate_m, only: readsulfate
60      use sugwd_m, only: sugwd      use sugwd_m, only: sugwd
61      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
62      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE temps, ONLY: annee_ref, day_ref, itau_phy
# Line 123  contains Line 124  contains
124      character(len = 6):: ocean = 'force '      character(len = 6):: ocean = 'force '
125      ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")      ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")
126    
     logical ok_ocean  
     SAVE ok_ocean  
   
127      ! "slab" ocean      ! "slab" ocean
128      REAL, save:: tslab(klon) ! temperature of ocean slab      REAL, save:: tslab(klon) ! temperature of ocean slab
129      REAL, save:: seaice(klon) ! glace de mer (kg/m2)      REAL, save:: seaice(klon) ! glace de mer (kg/m2)
# Line 169  contains Line 167  contains
167    
168      !MI Amip2 PV a theta constante      !MI Amip2 PV a theta constante
169    
170      INTEGER klevp1      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
171      PARAMETER(klevp1 = llm + 1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
   
     REAL swdn0(klon, klevp1), swdn(klon, klevp1)  
     REAL swup0(klon, klevp1), swup(klon, klevp1)  
172      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
173    
174      REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
175      REAL lwup0(klon, klevp1), lwup(klon, klevp1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
176      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
177    
178      !IM Amip2      !IM Amip2
# Line 270  contains Line 265  contains
265      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
266      ! soil temperature of surface fraction      ! soil temperature of surface fraction
267    
268      REAL fevap(klon, nbsrf)      REAL, save:: fevap(klon, nbsrf) ! evaporation
     SAVE fevap ! evaporation  
269      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
270      SAVE fluxlat      SAVE fluxlat
271    
# Line 320  contains Line 314  contains
314      SAVE qcondc      SAVE qcondc
315      REAL ema_work1(klon, llm), ema_work2(klon, llm)      REAL ema_work1(klon, llm), ema_work2(klon, llm)
316      SAVE ema_work1, ema_work2      SAVE ema_work1, ema_work2
317        REAL, save:: wd(klon)
     REAL wd(klon) ! sb  
     SAVE wd ! sb  
318    
319      ! Variables locales pour la couche limite (al1):      ! Variables locales pour la couche limite (al1):
320    
# Line 331  contains Line 323  contains
323      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
324      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
325    
326      !AA Pour phytrac      ! Pour phytrac :
327      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
328      REAL yu1(klon) ! vents dans la premiere couche U      REAL yu1(klon) ! vents dans la premiere couche U
329      REAL yv1(klon) ! vents dans la premiere couche V      REAL yv1(klon) ! vents dans la premiere couche V
# Line 355  contains Line 347  contains
347    
348      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
349    
350      REAL evap(klon), devap(klon) ! evaporation et sa derivee      REAL evap(klon), devap(klon) ! evaporation and its derivative
351      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
352      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
353      SAVE dlw      SAVE dlw
# Line 376  contains Line 368  contains
368      INTEGER julien      INTEGER julien
369    
370      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
371      REAL pctsrf(klon, nbsrf)      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
372      !IM      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE  
373    
     SAVE pctsrf ! sous-fraction du sol  
374      REAL albsol(klon)      REAL albsol(klon)
375      SAVE albsol ! albedo du sol total      SAVE albsol ! albedo du sol total
376      REAL albsollw(klon)      REAL albsollw(klon)
# Line 450  contains Line 440  contains
440      REAL dist, rmu0(klon), fract(klon)      REAL dist, rmu0(klon), fract(klon)
441      REAL zdtime ! pas de temps du rayonnement (s)      REAL zdtime ! pas de temps du rayonnement (s)
442      real zlongi      real zlongi
   
443      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
     LOGICAL zx_ajustq  
   
444      REAL za, zb      REAL za, zb
445      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zdelta, zcor
446      real zqsat(klon, llm)      real zqsat(klon, llm)
447      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
448      REAL t_coup      REAL, PARAMETER:: t_coup = 234.
     PARAMETER (t_coup = 234.0)  
   
449      REAL zphi(klon, llm)      REAL zphi(klon, llm)
450    
451      !IM cf. AM Variables locales pour la CLA (hbtm2)      !IM cf. AM Variables locales pour la CLA (hbtm2)
# Line 497  contains Line 482  contains
482      REAL rflag(klon) ! flag fonctionnement de convect      REAL rflag(klon) ! flag fonctionnement de convect
483      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
484      ! -- convect43:      ! -- convect43:
     INTEGER ntra ! nb traceurs pour convect4.3  
485      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)      REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
486      REAL dplcldt(klon), dplcldr(klon)      REAL dplcldt(klon), dplcldr(klon)
487    
# Line 515  contains Line 499  contains
499      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
500      REAL rneb(klon, llm)      REAL rneb(klon, llm)
501    
502      REAL pmfu(klon, llm), pmfd(klon, llm)      REAL mfu(klon, llm), mfd(klon, llm)
503      REAL pen_u(klon, llm), pen_d(klon, llm)      REAL pen_u(klon, llm), pen_d(klon, llm)
504      REAL pde_u(klon, llm), pde_d(klon, llm)      REAL pde_u(klon, llm), pde_d(klon, llm)
505      INTEGER kcbot(klon), kctop(klon), kdtop(klon)      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
# Line 579  contains Line 563  contains
563    
564      REAL zsto      REAL zsto
565    
     character(len = 20) modname  
     character(len = 80) abort_message  
566      logical ok_sync      logical ok_sync
567      real date0      real date0
568    
# Line 598  contains Line 580  contains
580      REAL ZRCPD      REAL ZRCPD
581    
582      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
583      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m      REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
584      REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille      REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
585      REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille      REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
586      !jq Aerosol effects (Johannes Quaas, 27/11/2003)  
587      REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]      ! Aerosol effects:
588    
589        REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
590    
591      REAL, save:: sulfate_pi(klon, llm)      REAL, save:: sulfate_pi(klon, llm)
592      ! (SO4 aerosol concentration, in ug/m3, pre-industrial value)      ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
593    
594      REAL cldtaupi(klon, llm)      REAL cldtaupi(klon, llm)
595      ! (Cloud optical thickness for pre-industrial (pi) aerosols)      ! cloud optical thickness for pre-industrial (pi) aerosols
596    
597      REAL re(klon, llm) ! Cloud droplet effective radius      REAL re(klon, llm) ! Cloud droplet effective radius
598      REAL fl(klon, llm) ! denominator of re      REAL fl(klon, llm) ! denominator of re
# Line 618  contains Line 602  contains
602      REAL, save:: cg_ae(klon, llm, 2)      REAL, save:: cg_ae(klon, llm, 2)
603    
604      REAL topswad(klon), solswad(klon) ! aerosol direct effect      REAL topswad(klon), solswad(klon) ! aerosol direct effect
     ! ok_ade --> ADE = topswad - topsw  
   
605      REAL topswai(klon), solswai(klon) ! aerosol indirect effect      REAL topswai(klon), solswai(klon) ! aerosol indirect effect
     ! ok_aie .and. ok_ade --> AIE = topswai - topswad  
     ! ok_aie .and. .not. ok_ade --> AIE = topswai - topsw  
606    
607      REAL aerindex(klon) ! POLDER aerosol index      REAL aerindex(klon) ! POLDER aerosol index
608    
# Line 630  contains Line 610  contains
610      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect      LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
611    
612      REAL:: bl95_b0 = 2., bl95_b1 = 0.2      REAL:: bl95_b0 = 2., bl95_b1 = 0.2
613      ! Parameters in the formula to link CDNC to aerosol mass conc      ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
614      ! (Boucher and Lohmann, 1995), used in nuage.F      ! B). They link cloud droplet number concentration to aerosol mass
615        ! concentration.
616    
617      SAVE u10m      SAVE u10m
618      SAVE v10m      SAVE v10m
# Line 662  contains Line 643  contains
643    
644      !----------------------------------------------------------------      !----------------------------------------------------------------
645    
646      modname = 'physiq'      IF (if_ebil >= 1) zero_v = 0.
     IF (if_ebil >= 1) THEN  
        DO i = 1, klon  
           zero_v(i) = 0.  
        END DO  
     END IF  
647      ok_sync = .TRUE.      ok_sync = .TRUE.
648      IF (nqmx < 2) THEN      IF (nqmx < 2) CALL abort_gcm('physiq', &
649         abort_message = 'eaux vapeur et liquide sont indispensables'           'eaux vapeur et liquide sont indispensables', 1)
        CALL abort_gcm(modname, abort_message, 1)  
     ENDIF  
650    
651      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
652         ! initialiser         ! initialiser
# Line 687  contains Line 661  contains
661         cg_ae = 0.         cg_ae = 0.
662         rain_con(:) = 0.         rain_con(:) = 0.
663         snow_con(:) = 0.         snow_con(:) = 0.
        bl95_b0 = 0.  
        bl95_b1 = 0.  
664         topswai(:) = 0.         topswai(:) = 0.
665         topswad(:) = 0.         topswad(:) = 0.
666         solswai(:) = 0.         solswai(:) = 0.
# Line 720  contains Line 692  contains
692         read(unit=*, nml=physiq_nml)         read(unit=*, nml=physiq_nml)
693         write(unit_nml, nml=physiq_nml)         write(unit_nml, nml=physiq_nml)
694    
        ! Appel à la lecture du run.def physique  
695         call conf_phys         call conf_phys
696    
697         ! Initialiser les compteurs:         ! Initialiser les compteurs:
# Line 735  contains Line 706  contains
706              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0)
707    
708         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
709         q2 = 1.e-8         q2 = 1e-8
710    
711         radpas = NINT(86400. / dtphys / nbapp_rad)         radpas = NINT(86400. / dtphys / nbapp_rad)
712    
# Line 743  contains Line 714  contains
714         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
715    
716         PRINT *, 'cycle_diurne = ', cycle_diurne         PRINT *, 'cycle_diurne = ', cycle_diurne
717           CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &
718                ok_instan, ok_region)
719    
720         IF(ocean.NE.'force ') THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
           ok_ocean = .TRUE.  
        ENDIF  
   
        CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &  
             ok_region)  
   
        IF (dtphys*REAL(radpas) > 21600..AND.cycle_diurne) THEN  
           print *, 'Nbre d appels au rayonnement insuffisant'  
721            print *, "Au minimum 4 appels par jour si cycle diurne"            print *, "Au minimum 4 appels par jour si cycle diurne"
722            abort_message = 'Nbre d appels au rayonnement insuffisant'            call abort_gcm('physiq', &
723            call abort_gcm(modname, abort_message, 1)                 "Nombre d'appels au rayonnement insuffisant", 1)
724         ENDIF         ENDIF
        print *, "Clef pour la convection, iflag_con = ", iflag_con  
725    
726         ! Initialisation pour la convection de K.E. (sb):         ! Initialisation pour le schéma de convection d'Emanuel :
727         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
728            print *, "Convection de Kerry Emanuel 4.3"            ibas_con = 1
729              itop_con = 1
           DO i = 1, klon  
              ibas_con(i) = 1  
              itop_con(i) = 1  
           ENDDO  
730         ENDIF         ENDIF
731    
732         IF (ok_orodr) THEN         IF (ok_orodr) THEN
# Line 796  contains Line 756  contains
756         call ini_histday(dtphys, ok_journe, nid_day, nqmx)         call ini_histday(dtphys, ok_journe, nid_day, nqmx)
757         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
758         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
759         !XXXPB Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
760         WRITE(*, *) 'physiq date0: ', date0         print *, 'physiq date0: ', date0
761      ENDIF test_firstcal      ENDIF test_firstcal
762    
763      ! Mettre a zero des variables de sortie (pour securite)      ! Mettre a zero des variables de sortie (pour securite)
764    
765      DO i = 1, klon      DO i = 1, klon
766         d_ps(i) = 0.0         d_ps(i) = 0.
767      ENDDO      ENDDO
768      DO iq = 1, nqmx      DO iq = 1, nqmx
769         DO k = 1, llm         DO k = 1, llm
770            DO i = 1, klon            DO i = 1, klon
771               d_qx(i, k, iq) = 0.0               d_qx(i, k, iq) = 0.
772            ENDDO            ENDDO
773         ENDDO         ENDDO
774      ENDDO      ENDDO
# Line 893  contains Line 853  contains
853    
854      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
855    
856      ! Mettre en action les conditions aux limites (albedo, sst, etc.).      ! Mettre en action les conditions aux limites (albedo, sst etc.).
857    
858      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
859      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
# Line 962  contains Line 922  contains
922      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
923         DO i = 1, klon         DO i = 1, klon
924            fsollw(i, nsrf) = sollw(i) &            fsollw(i, nsrf) = sollw(i) &
925                 + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))
926            fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))
927         ENDDO         ENDDO
928      ENDDO      ENDDO
929    
# Line 991  contains Line 951  contains
951      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
952         DO k = 1, llm         DO k = 1, llm
953            DO i = 1, klon            DO i = 1, klon
954               zxfluxt(i, k) = zxfluxt(i, k) + &               zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
955                    fluxt(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
956               zxfluxq(i, k) = zxfluxq(i, k) + &               zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
957                    fluxq(i, k, nsrf) * pctsrf(i, nsrf)               zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
              zxfluxu(i, k) = zxfluxu(i, k) + &  
                   fluxu(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxv(i, k) = zxfluxv(i, k) + &  
                   fluxv(i, k, nsrf) * pctsrf(i, nsrf)  
958            END DO            END DO
959         END DO         END DO
960      END DO      END DO
961      DO i = 1, klon      DO i = 1, klon
962         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
963         evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'évaporation au sol
964         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
965      ENDDO      ENDDO
966    
# Line 1051  contains Line 1007  contains
1007         s_trmb2(i) = 0.0         s_trmb2(i) = 0.0
1008         s_trmb3(i) = 0.0         s_trmb3(i) = 0.0
1009    
1010         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
1011              pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.)  >  EPSFRA) &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
1012              THEN              'physiq : problème sous surface au point ', i, pctsrf(i, 1 : nbsrf)
           WRITE(*, *) 'physiq : pb sous surface au point ', i, &  
                pctsrf(i, 1 : nbsrf)  
        ENDIF  
1013      ENDDO      ENDDO
1014      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1015         DO i = 1, klon         DO i = 1, klon
# Line 1113  contains Line 1066  contains
1066      ! Calculer la derive du flux infrarouge      ! Calculer la derive du flux infrarouge
1067    
1068      DO i = 1, klon      DO i = 1, klon
1069         dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
1070      ENDDO      ENDDO
1071    
1072      ! Appeler la convection (au choix)      ! Appeler la convection (au choix)
1073    
1074      DO k = 1, llm      DO k = 1, llm
1075         DO i = 1, klon         DO i = 1, klon
1076            conv_q(i, k) = d_q_dyn(i, k) &            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys
1077                 + d_q_vdf(i, k)/dtphys            conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys
           conv_t(i, k) = d_t_dyn(i, k) &  
                + d_t_vdf(i, k)/dtphys  
1078         ENDDO         ENDDO
1079      ENDDO      ENDDO
1080    
1081      IF (check) THEN      IF (check) THEN
1082         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1083         print *, "avantcon = ", za         print *, "avantcon = ", za
1084      ENDIF      ENDIF
     zx_ajustq = iflag_con == 2  
     IF (zx_ajustq) THEN  
        DO i = 1, klon  
           z_avant(i) = 0.0  
        ENDDO  
        DO k = 1, llm  
           DO i = 1, klon  
              z_avant(i) = z_avant(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
     ENDIF  
1085    
1086      select case (iflag_con)      if (iflag_con == 2) then
1087      case (2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
1088         CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
1089              zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, &              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
1090              pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, &              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
1091              pmflxs)              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
1092                kdtop, pmflxr, pmflxs)
1093         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
1094         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
1095         DO i = 1, klon         ibas_con = llm + 1 - kcbot
1096            ibas_con(i) = llm + 1 - kcbot(i)         itop_con = llm + 1 - kctop
1097            itop_con(i) = llm + 1 - kctop(i)      else
1098         ENDDO         ! iflag_con >= 3
1099      case (3:)         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &
1100         ! number of tracers for the convection scheme of Kerry Emanuel:              v_seri, tr_seri, ema_work1, ema_work2, d_t_con, d_q_con, &
1101                d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &
1102                itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &
1103                pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &
1104                wd, pmflxr, pmflxs, da, phi, mp, ntra=1)
1105           ! (number of tracers for the convection scheme of Kerry Emanuel:
1106         ! la partie traceurs est faite dans phytrac         ! la partie traceurs est faite dans phytrac
1107         ! on met ntra = 1 pour limiter les appels mais on peut         ! on met ntra = 1 pour limiter les appels mais on peut
1108         ! supprimer les calculs / ftra.         ! supprimer les calculs / ftra.)
        ntra = 1  
        ! Schéma de convection modularisé et vectorisé :  
        ! (driver commun aux versions 3 et 4)  
   
        CALL concvl(iflag_con, dtphys, paprs, play, t_seri, q_seri, u_seri, &  
             v_seri, tr_seri, ntra, ema_work1, ema_work2, d_t_con, d_q_con, &  
             d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, itop_con, &  
             upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, pbase, bbase, &  
             dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, pmflxr, pmflxs, &  
             da, phi, mp)  
        clwcon0 = qcondc  
        pmfu = upwd + dnwd  
1109    
1110         IF (.NOT. ok_gust) THEN         clwcon0 = qcondc
1111            do i = 1, klon         mfu = upwd + dnwd
1112               wd(i) = 0.0         IF (.NOT. ok_gust) wd = 0.
           enddo  
        ENDIF  
1113    
1114         ! Calcul des propriétés des nuages convectifs         ! Calcul des propriétés des nuages convectifs
1115    
# Line 1186  contains Line 1118  contains
1118               zx_t = t_seri(i, k)               zx_t = t_seri(i, k)
1119               IF (thermcep) THEN               IF (thermcep) THEN
1120                  zdelta = MAX(0., SIGN(1., rtt-zx_t))                  zdelta = MAX(0., SIGN(1., rtt-zx_t))
1121                  zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)                  zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)
1122                  zx_qs = MIN(0.5, zx_qs)                  zx_qs = MIN(0.5, zx_qs)
1123                  zcor = 1./(1.-retv*zx_qs)                  zcor = 1./(1.-retv*zx_qs)
1124                  zx_qs = zx_qs*zcor                  zx_qs = zx_qs*zcor
# Line 1202  contains Line 1134  contains
1134         ENDDO         ENDDO
1135    
1136         ! calcul des proprietes des nuages convectifs         ! calcul des proprietes des nuages convectifs
1137         clwcon0 = fact_cldcon*clwcon0         clwcon0 = fact_cldcon * clwcon0
1138         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
1139              rnebcon0)              rnebcon0)
1140      case default      END if
        print *, "iflag_con non-prevu", iflag_con  
        stop 1  
     END select  
1141    
1142      DO k = 1, llm      DO k = 1, llm
1143         DO i = 1, klon         DO i = 1, klon
# Line 1242  contains Line 1171  contains
1171         zx_t = zx_t/za*dtphys         zx_t = zx_t/za*dtphys
1172         print *, "Precip = ", zx_t         print *, "Precip = ", zx_t
1173      ENDIF      ENDIF
1174      IF (zx_ajustq) THEN  
1175         DO i = 1, klon      IF (iflag_con == 2) THEN
1176            z_apres(i) = 0.0         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1177         ENDDO         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
        DO k = 1, llm  
           DO i = 1, klon  
              z_apres(i) = z_apres(i) + (q_seri(i, k) + ql_seri(i, k)) &  
                   *zmasse(i, k)  
           ENDDO  
        ENDDO  
        DO i = 1, klon  
           z_factor(i) = (z_avant(i)-(rain_con(i) + snow_con(i))*dtphys) &  
                /z_apres(i)  
        ENDDO  
1178         DO k = 1, llm         DO k = 1, llm
1179            DO i = 1, klon            DO i = 1, klon
1180               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN               IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
# Line 1264  contains Line 1183  contains
1183            ENDDO            ENDDO
1184         ENDDO         ENDDO
1185      ENDIF      ENDIF
     zx_ajustq = .FALSE.  
1186    
1187      ! Convection sèche (thermiques ou ajustement)      ! Convection sèche (thermiques ou ajustement)
1188    
# Line 1295  contains Line 1213  contains
1213    
1214      ! Caclul des ratqs      ! Caclul des ratqs
1215    
1216      ! ratqs convectifs a l'ancienne en fonction de q(z = 0)-q / q      ! ratqs convectifs à l'ancienne en fonction de (q(z = 0) - q) / q
1217      ! on ecrase le tableau ratqsc calcule par clouds_gno      ! on écrase le tableau ratqsc calculé par clouds_gno
1218      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
1219         do k = 1, llm         do k = 1, llm
1220            do i = 1, klon            do i = 1, klon
1221               if(ptconv(i, k)) then               if(ptconv(i, k)) then
1222                  ratqsc(i, k) = ratqsbas &                  ratqsc(i, k) = ratqsbas + fact_cldcon &
1223                       +fact_cldcon*(q_seri(i, 1)-q_seri(i, k))/q_seri(i, k)                       * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
1224               else               else
1225                  ratqsc(i, k) = 0.                  ratqsc(i, k) = 0.
1226               endif               endif
# Line 1313  contains Line 1231  contains
1231      ! ratqs stables      ! ratqs stables
1232      do k = 1, llm      do k = 1, llm
1233         do i = 1, klon         do i = 1, klon
1234            ratqss(i, k) = ratqsbas + (ratqshaut-ratqsbas)* &            ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1235                 min((paprs(i, 1)-play(i, k))/(paprs(i, 1)-30000.), 1.)                 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1236         enddo         enddo
1237      enddo      enddo
1238    
1239      ! ratqs final      ! ratqs final
1240      if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then      if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1241         ! les ratqs sont une conbinaison de ratqss et ratqsc         ! les ratqs sont une conbinaison de ratqss et ratqsc
1242         ! ratqs final         ! ratqs final
1243         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de         ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1244         ! relaxation des ratqs         ! relaxation des ratqs
1245         facteur = exp(-dtphys*facttemps)         ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
        ratqs = max(ratqs*facteur, ratqss)  
1246         ratqs = max(ratqs, ratqsc)         ratqs = max(ratqs, ratqsc)
1247      else      else
1248         ! on ne prend que le ratqs stable pour fisrtilp         ! on ne prend que le ratqs stable pour fisrtilp
# Line 1396  contains Line 1313  contains
1313         endif         endif
1314    
1315         ! Nuages diagnostiques pour Tiedtke         ! Nuages diagnostiques pour Tiedtke
1316         CALL diagcld1(paprs, play, &         CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1317              rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &              itop_con, diafra, dialiq)
             diafra, dialiq)  
1318         DO k = 1, llm         DO k = 1, llm
1319            DO i = 1, klon            DO i = 1, klon
1320               IF (diafra(i, k) > cldfra(i, k)) THEN               IF (diafra(i, k) > cldfra(i, k)) THEN
# Line 1414  contains Line 1330  contains
1330         facteur = dtphys *facttemps         facteur = dtphys *facttemps
1331         do k = 1, llm         do k = 1, llm
1332            do i = 1, klon            do i = 1, klon
1333               rnebcon(i, k) = rnebcon(i, k)*facteur               rnebcon(i, k) = rnebcon(i, k) * facteur
1334               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &               if (rnebcon0(i, k)*clwcon0(i, k) > rnebcon(i, k)*clwcon(i, k)) &
1335                    then                    then
1336                  rnebcon(i, k) = rnebcon0(i, k)                  rnebcon(i, k) = rnebcon0(i, k)
# Line 1490  contains Line 1406  contains
1406    
1407      ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :      ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :
1408      if (ok_newmicro) then      if (ok_newmicro) then
1409         CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1410              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1411              fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, &              sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
             re, fl)  
1412      else      else
1413         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1414              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
# Line 1637  contains Line 1552  contains
1552    
1553      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1554      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &
1555           dtphys, u, t, paprs, play, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &           dtphys, u, t, paprs, play, mfu, mfd, pen_u, pde_u, pen_d, pde_d, &
1556           ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &           ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &
1557           frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &           frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &
1558           pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)           pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)
1559    
1560      IF (offline) THEN      IF (offline) THEN
1561         call phystokenc(dtphys, rlon, rlat, t, pmfu, pmfd, pen_u, pde_u, &         call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, pde_u, &
1562              pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &              pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1563              pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)              pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
1564      ENDIF      ENDIF
# Line 1689  contains Line 1604  contains
1604    
1605      ! SORTIES      ! SORTIES
1606    
1607      !cc prw = eau precipitable      ! prw = eau precipitable
1608      DO i = 1, klon      DO i = 1, klon
1609         prw(i) = 0.         prw(i) = 0.
1610         DO k = 1, llm         DO k = 1, llm
Legend:
Removed from v.68  
changed lines
  Added in v.71
  ViewVC Help
Powered by ViewVC 1.1.21