/[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 97 by guez, Fri Apr 25 14:58:31 2014 UTC revision 99 by guez, Wed Jul 2 18:39:15 2014 UTC
# Line 35  contains Line 35  contains
35      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
36      use diagphy_m, only: diagphy      use diagphy_m, only: diagphy
37      USE dimens_m, ONLY: llm, nqmx      USE dimens_m, ONLY: llm, nqmx
38      USE dimphy, ONLY: klon, nbtr      USE dimphy, ONLY: klon
39      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
40      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
41      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
# Line 45  contains Line 45  contains
45           nbsrf           nbsrf
46      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins
47      use newmicro_m, only: newmicro      use newmicro_m, only: newmicro
48      USE oasis_m, ONLY: ok_oasis      USE orbite_m, ONLY: orbite
     USE orbite_m, ONLY: orbite, zenang  
49      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
50      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
51      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
# Line 61  contains Line 60  contains
60      use unit_nml_m, only: unit_nml      use unit_nml_m, only: unit_nml
61      USE ymds2ju_m, ONLY: ymds2ju      USE ymds2ju_m, ONLY: ymds2ju
62      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
63        use zenang_m, only: zenang
64    
65      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
66    
# Line 70  contains Line 70  contains
70      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
71      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)
72    
73      REAL, intent(in):: paprs(klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
74      ! (pression pour chaque inter-couche, en Pa)      ! pression pour chaque inter-couche, en Pa
75    
76      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
77      ! (input pression pour le mileu de chaque couche (en Pa))      ! pression pour le mileu de chaque couche (en Pa)
78    
79      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
80      ! géopotentiel de chaque couche (référence sol)      ! géopotentiel de chaque couche (référence sol)
81    
82      REAL, intent(in):: pphis(klon) ! géopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
83    
84      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
85      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m/s
86    
87      REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s
88      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
89    
90      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
91      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
92    
93      REAL, intent(in):: omega(klon, llm) ! vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
94      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m s-2)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
95      REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m s-2)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
96      REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)
97      REAL, intent(out):: d_qx(klon, llm, nqmx) ! tendance physique de "qx" (s-1)  
98        REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
99        ! tendance physique de "qx" (s-1)
100    
101      ! Local:      ! Local:
102    
103      LOGICAL:: firstcal = .true.      LOGICAL:: firstcal = .true.
104    
     INTEGER nbteta  
     PARAMETER(nbteta = 3)  
   
105      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
106      PARAMETER (ok_gust = .FALSE.)      PARAMETER (ok_gust = .FALSE.)
107    
108      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL, PARAMETER:: check = .FALSE.
109      PARAMETER (check = .FALSE.)      ! Verifier la conservation du modele en eau
110    
111      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
112      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
113    
     ! Parametres lies au coupleur OASIS:  
     INTEGER, SAVE:: npas, nexca  
     logical rnpb  
     parameter(rnpb = .true.)  
   
     character(len = 6):: ocean = 'force '  
     ! (type de mod\`ele oc\'ean \`a utiliser: "force" ou "slab" mais  
     ! pas "couple")  
   
114      ! "slab" ocean      ! "slab" ocean
115      REAL, save:: tslab(klon) ! temperature of ocean slab      REAL, save:: tslab(klon) ! temperature of ocean slab
116      REAL, save:: seaice(klon) ! glace de mer (kg/m2)      REAL, save:: seaice(klon) ! glace de mer (kg/m2)
117      REAL fluxo(klon) ! flux turbulents ocean-glace de mer      REAL fluxo(klon) ! flux turbulents ocean-glace de mer
118      REAL fluxg(klon) ! flux turbulents ocean-atmosphere      REAL fluxg(klon) ! flux turbulents ocean-atmosphere
119    
     ! Modele thermique du sol, a activer pour le cycle diurne:  
     logical:: ok_veget = .false. ! type de modele de vegetation utilise  
   
120      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
121      ! sorties journalieres, mensuelles et instantanees dans les      ! sorties journalieres, mensuelles et instantanees dans les
122      ! fichiers histday, histmth et histins      ! fichiers histday, histmth et histins
# Line 142  contains Line 129  contains
129      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
130      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
131    
132      INTEGER ivap ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
133      PARAMETER (ivap = 1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     INTEGER iliq ! indice de traceurs pour eau liquide  
     PARAMETER (iliq = 2)  
134    
135      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
136      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
# Line 155  contains Line 140  contains
140    
141      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
142    
     ! Amip2 PV a theta constante  
   
     CHARACTER(LEN = 3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
   
     ! Amip2 PV a theta constante  
   
143      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
144      REAL swup0(klon, llm + 1), swup(klon, llm + 1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
145      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
# Line 263  contains Line 239  contains
239      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
240      SAVE fluxlat      SAVE fluxlat
241    
242      REAL fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
243      SAVE fqsurf ! humidite de l'air au contact de la surface      ! humidite de l'air au contact de la surface
244    
245      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol
246        REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
247      REAL fsnow(klon, nbsrf)      REAL, save:: falbe(klon, nbsrf) ! albedo par type de surface
248      SAVE fsnow ! epaisseur neigeuse      REAL, save:: falblw(klon, nbsrf) ! albedo par type de surface
   
     REAL falbe(klon, nbsrf)  
     SAVE falbe ! albedo par type de surface  
     REAL falblw(klon, nbsrf)  
     SAVE falblw ! albedo par type de surface  
249    
250      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
251      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
# Line 348  contains Line 319  contains
319      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
320      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
321    
322      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
323      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
324    
325      ! Conditions aux limites      ! Conditions aux limites
326    
327      INTEGER julien      INTEGER julien
   
328      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
329      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
330      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
331        REAL, save:: albsol(klon) ! albedo du sol total
332      REAL albsol(klon)      REAL, save:: albsollw(klon) ! albedo du sol total
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
333      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
334    
335      ! Declaration des procedures appelees      ! Declaration des procedures appelees
# Line 423  contains Line 388  contains
388      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
389    
390      REAL dist, rmu0(klon), fract(klon)      REAL dist, rmu0(klon), fract(klon)
     REAL zdtime ! pas de temps du rayonnement (s)  
391      real zlongi      real zlongi
392      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
393      REAL za, zb      REAL za, zb
# Line 512  contains Line 476  contains
476      ! Variables locales pour effectuer les appels en s\'erie :      ! Variables locales pour effectuer les appels en s\'erie :
477    
478      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
479      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
480      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
481      REAL tr_seri(klon, llm, nbtr)      REAL tr_seri(klon, llm, nqmx - 2)
482    
483      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
484    
# Line 537  contains Line 501  contains
501    
502      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
503      REAL ztsol(klon)      REAL ztsol(klon)
504      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
505      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
506      REAL zero_v(klon)      REAL zero_v(klon)
507      CHARACTER(LEN = 20) tit      CHARACTER(LEN = 20) tit
508      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
# Line 603  contains Line 566  contains
566    
567      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
568    
569      namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &      namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
570           fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &           facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
571           ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &           ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals
          nsplit_thermals  
572    
573      !----------------------------------------------------------------      !----------------------------------------------------------------
574    
# Line 625  contains Line 587  contains
587         piz_ae = 0.         piz_ae = 0.
588         tau_ae = 0.         tau_ae = 0.
589         cg_ae = 0.         cg_ae = 0.
590         rain_con(:) = 0.         rain_con = 0.
591         snow_con(:) = 0.         snow_con = 0.
592         topswai(:) = 0.         topswai = 0.
593         topswad(:) = 0.         topswad = 0.
594         solswai(:) = 0.         solswai = 0.
595         solswad(:) = 0.         solswad = 0.
596    
597         d_u_con = 0.         d_u_con = 0.
598         d_v_con = 0.         d_v_con = 0.
# Line 665  contains Line 627  contains
627         frugs = 0.         frugs = 0.
628         itap = 0         itap = 0
629         itaprad = 0         itaprad = 0
630         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &         CALL phyetat0(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &
631              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              fsnow, falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollw, &
632              snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &              dlw, radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, &
633              zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &              zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
634              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)              run_off_lic_0, sig1, w01)
635    
636         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
637         q2 = 1e-8         q2 = 1e-8
# Line 680  contains Line 642  contains
642         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
643    
644         PRINT *, 'cycle_diurne = ', cycle_diurne         PRINT *, 'cycle_diurne = ', cycle_diurne
645         CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &         CALL printflag(radpas, ok_journe, ok_instan, ok_region)
             ok_instan, ok_region)  
646    
647         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
648            print *, "Au minimum 4 appels par jour si cycle diurne"            print *, "Au minimum 4 appels par jour si cycle diurne"
# Line 711  contains Line 672  contains
672         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
673         ecrit_reg = NINT(ecrit_reg/dtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
674    
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
675         ! Initialisation des sorties         ! Initialisation des sorties
676    
677         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
# Line 724  contains Line 680  contains
680         print *, 'physiq date0: ', date0         print *, 'physiq date0: ', date0
681      ENDIF test_firstcal      ENDIF test_firstcal
682    
     ! Mettre a zero des variables de sortie (pour securite)  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
683      ! We will modify variables *_seri and we will not touch variables      ! We will modify variables *_seri and we will not touch variables
684      ! u, v, h, q:      ! u, v, t, qx:
685      DO k = 1, llm      t_seri = t
686         DO i = 1, klon      u_seri = u
687            t_seri(i, k) = t(i, k)      v_seri = v
688            u_seri(i, k) = u(i, k)      q_seri = qx(:, :, ivap)
689            v_seri(i, k) = v(i, k)      ql_seri = qx(:, :, iliq)
690            q_seri(i, k) = qx(i, k, ivap)      tr_seri = qx(:, :, 3: nqmx)
           ql_seri(i, k) = qx(i, k, iliq)  
           qs_seri(i, k) = 0.  
        ENDDO  
     ENDDO  
     IF (nqmx >= 3) THEN  
        tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)  
     ELSE  
        tr_seri(:, :, 1) = 0.  
     ENDIF  
691    
692      DO i = 1, klon      ztsol = sum(ftsol * pctsrf, dim = 2)
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
693    
694      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
695         tit = 'after dynamics'         tit = 'after dynamics'
696         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
697              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
698         ! Comme les tendances de la physique sont ajout\'es dans la         ! Comme les tendances de la physique sont ajout\'es dans la
699         !  dynamique, la variation d'enthalpie par la dynamique devrait         !  dynamique, la variation d'enthalpie par la dynamique devrait
700         !  \^etre \'egale \`a la variation de la physique au pas de temps         !  \^etre \'egale \`a la variation de la physique au pas de temps
# Line 768  contains Line 702  contains
702         !  nulle.         !  nulle.
703         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
704              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
705              d_qt, 0., fs_bound, fq_bound)              d_qt, 0.)
706      END IF      END IF
707    
708      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
# Line 799  contains Line 733  contains
733      ! Check temperatures:      ! Check temperatures:
734      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
735    
736      ! Incrementer le compteur de la physique      ! Incrémenter le compteur de la physique
737      itap = itap + 1      itap = itap + 1
738      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(NINT(rdayvrai), 360)
739      if (julien == 0) julien = 360      if (julien == 0) julien = 360
740    
741      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
742    
743      ! Mettre en action les conditions aux limites (albedo, sst etc.).      ! Prescrire l'ozone :
   
     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.  
744      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
745    
746      ! \'Evaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
# Line 825  contains Line 757  contains
757      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
758         tit = 'after reevap'         tit = 'after reevap'
759         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
760              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
761         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
762              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
   
763      END IF      END IF
764    
765      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
766        zxrugs = sum(frugs * pctsrf, dim = 2)
     DO i = 1, klon  
        zxrugs(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)  
        ENDDO  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
767    
768      ! calculs necessaires au calcul de l'albedo dans l'interface      ! Calculs nécessaires au calcul de l'albedo dans l'interface
769    
770      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), zlongi, dist)
771      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
772         zdtime = dtphys * REAL(radpas)         CALL zenang(zlongi, time, dtphys * REAL(radpas), rmu0, fract)
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
773      ELSE      ELSE
774         rmu0 = -999.999         rmu0 = -999.999
775      ENDIF      ENDIF
776    
777      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
778      albsol(:) = 0.      albsol = sum(falbe * pctsrf, dim = 2)
779      albsollw(:) = 0.      albsollw = sum(falblw * pctsrf, dim = 2)
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)  
           albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
780    
781      ! R\'epartition sous maille des flux longwave et shortwave      ! R\'epartition sous maille des flux longwave et shortwave
782      ! R\'epartition du longwave par sous-surface lin\'earis\'ee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
783    
784      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
785         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
786            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
787                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
788            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))      END forall
        ENDDO  
     ENDDO  
789    
790      fder = dlw      fder = dlw
791    
792      ! Couche limite:      ! Couche limite:
793    
794      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
795           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &           v_seri, julien, rmu0, co2_ppm, ftsol, soil_model, &
796           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, &
797           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           fsnow, fqsurf, fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, &
798           rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &           fsolsw, fsollw, fder, rlat, frugs, firstcal, agesno, rugoro, &
799           frugs, firstcal, agesno, rugoro, d_t_vdf, &           d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, &
800           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, &
801           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, &
802           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           trmb3, plcl, fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, &
803           fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)           seaice)
804    
805      ! Incr\'ementation des flux      ! Incr\'ementation des flux
806    
# Line 929  contains Line 836  contains
836      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
837         tit = 'after clmain'         tit = 'after clmain'
838         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
839              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
840         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
841              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
842      END IF      END IF
843    
844      ! Update surface temperature:      ! Update surface temperature:
# Line 1016  contains Line 921  contains
921         ENDDO         ENDDO
922      ENDDO      ENDDO
923    
924      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
925    
926      DO i = 1, klon      DO i = 1, klon
927         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
# Line 1026  contains Line 931  contains
931    
932      DO k = 1, llm      DO k = 1, llm
933         DO i = 1, klon         DO i = 1, klon
934            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys            conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k) / dtphys
935            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
936         ENDDO         ENDDO
937      ENDDO      ENDDO
938    
939      IF (check) THEN      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
940    
941      if (iflag_con == 2) then      if (iflag_con == 2) then
942         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
# Line 1050  contains Line 952  contains
952      else      else
953         ! iflag_con >= 3         ! iflag_con >= 3
954    
955           da = 0.
956           mp = 0.
957           phi = 0.
958         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
959              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &
960              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &
# Line 1101  contains Line 1006  contains
1006      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1007         tit = 'after convect'         tit = 'after convect'
1008         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1009              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1010         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1011              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
1012      END IF      END IF
1013    
1014      IF (check) THEN      IF (check) THEN
1015         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1016         print *, "aprescon = ", za         print *, "aprescon = ", za
1017         zx_t = 0.         zx_t = 0.
1018         za = 0.         za = 0.
# Line 1157  contains Line 1060  contains
1060      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1061         tit = 'after dry_adjust'         tit = 'after dry_adjust'
1062         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1063              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1064      END IF      END IF
1065    
1066      ! Caclul des ratqs      ! Caclul des ratqs
# Line 1216  contains Line 1118  contains
1118         ENDDO         ENDDO
1119      ENDDO      ENDDO
1120      IF (check) THEN      IF (check) THEN
1121         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1122         print *, "apresilp = ", za         print *, "apresilp = ", za
1123         zx_t = 0.         zx_t = 0.
1124         za = 0.         za = 0.
# Line 1232  contains Line 1134  contains
1134      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1135         tit = 'after fisrt'         tit = 'after fisrt'
1136         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1137              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1138         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1139              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
1140      END IF      END IF
1141    
1142      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
# Line 1313  contains Line 1213  contains
1213      ENDDO      ENDDO
1214    
1215      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1216           dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1217           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1218    
1219      ! Humidit\'e relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
1220      DO k = 1, llm      DO k = 1, llm
# Line 1398  contains Line 1298  contains
1298      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1299         tit = 'after rad'         tit = 'after rad'
1300         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1301              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1302         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1303              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
1304      END IF      END IF
1305    
1306      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
# Line 1496  contains Line 1394  contains
1394           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1395    
1396      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1397           2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1398           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1399    
1400      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1401      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, u, t, &
1402           dtphys, u, t, paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, &           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1403           entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, albsol, rhcl, &
1404           albsol, rhcl, cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, &           cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, &
1405           psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)           mp, upwd, dnwd, tr_seri, zmasse)
1406    
1407      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1408           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
# Line 1535  contains Line 1433  contains
1433      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1434         tit = 'after physic'         tit = 'after physic'
1435         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1436              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1437         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1438         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1439         ! est egale a la variation de la physique au pas de temps precedent.         ! est egale a la variation de la physique au pas de temps precedent.
1440         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1441         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1442              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
   
1443         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1444      END IF      END IF
1445    
1446      ! SORTIES      ! SORTIES
# Line 1571  contains Line 1465  contains
1465         ENDDO         ENDDO
1466      ENDDO      ENDDO
1467    
1468      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1469         DO iq = 3, nqmx         DO k = 1, llm
1470            DO k = 1, llm            DO i = 1, klon
1471               DO i = 1, klon               d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys
                 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys  
              ENDDO  
1472            ENDDO            ENDDO
1473         ENDDO         ENDDO
1474      ENDIF      ENDDO
1475    
1476      ! Sauvegarder les valeurs de t et q a la fin de la physique:      ! Sauvegarder les valeurs de t et q a la fin de la physique:
1477      DO k = 1, llm      DO k = 1, llm
Legend:
Removed from v.97  
changed lines
  Added in v.99
  ViewVC Help
Powered by ViewVC 1.1.21