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

Diff of /trunk/phylmd/physiq.f

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

trunk/phylmd/physiq.f revision 101 by guez, Mon Jul 7 17:45:21 2014 UTC trunk/Sources/phylmd/physiq.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC
# Line 4  module physiq_m Line 4  module physiq_m
4    
5  contains  contains
6    
7    SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &    SUBROUTINE physiq(lafin, dayvrai, time, dtphys, paprs, play, pphi, pphis, &
8         u, v, t, qx, omega, d_u, d_v, d_t, d_qx)         u, v, t, qx, omega, d_u, d_v, d_t, d_qx)
9    
10      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
# Line 38  contains Line 38  contains
38      USE dimphy, ONLY: klon      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 dynetat0_m, only: day_ref, annee_ref
42      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
43      use fisrtilp_m, only: fisrtilp      use fisrtilp_m, only: fisrtilp
44      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
# Line 54  contains Line 55  contains
55      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
56      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
57      use readsulfate_m, only: readsulfate      use readsulfate_m, only: readsulfate
58        use readsulfate_preind_m, only: readsulfate_preind
59      use sugwd_m, only: sugwd      use sugwd_m, only: sugwd
60      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
61      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE temps, ONLY: itau_phy
62      use unit_nml_m, only: unit_nml      use unit_nml_m, only: unit_nml
63      USE ymds2ju_m, ONLY: ymds2ju      USE ymds2ju_m, ONLY: ymds2ju
64      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
# Line 64  contains Line 66  contains
66    
67      logical, intent(in):: lafin ! dernier passage      logical, intent(in):: lafin ! dernier passage
68    
69      REAL, intent(in):: rdayvrai      integer, intent(in):: dayvrai
70      ! (elapsed time since January 1st 0h of the starting year, in days)      ! current day number, based at value 1 on January 1st of annee_ref
71    
72      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
73      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)      REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)
# Line 222  contains Line 224  contains
224      ! Variables propres a la physique      ! Variables propres a la physique
225    
226      INTEGER, save:: radpas      INTEGER, save:: radpas
227      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
228      ! "physiq".)      ! "physiq".
229    
230      REAL radsol(klon)      REAL radsol(klon)
231      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
# Line 392  contains Line 394  contains
394    
395      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
396    
397      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
398      real zlongi      real longi
399      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
400      REAL za, zb      REAL za, zb
401      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
402      real zqsat(klon, llm)      real zqsat(klon, llm)
403      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
404      REAL, PARAMETER:: t_coup = 234.      REAL, PARAMETER:: t_coup = 234.
# Line 646  contains Line 648  contains
648         ! on remet le calendrier a zero         ! on remet le calendrier a zero
649         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
650    
        PRINT *, 'cycle_diurne = ', cycle_diurne  
651         CALL printflag(radpas, ok_journe, ok_instan, ok_region)         CALL printflag(radpas, ok_journe, ok_instan, ok_region)
652    
653         IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN         IF (dtphys * radpas > 21600. .AND. cycle_diurne) THEN
654            print *, "Au minimum 4 appels par jour si cycle diurne"            print *, "Au minimum 4 appels par jour si cycle diurne"
655            call abort_gcm('physiq', &            call abort_gcm('physiq', &
656                 "Nombre d'appels au rayonnement insuffisant", 1)                 "Nombre d'appels au rayonnement insuffisant", 1)
# Line 680  contains Line 681  contains
681         ! Initialisation des sorties         ! Initialisation des sorties
682    
683         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
684         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
685         ! Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
686         print *, 'physiq date0: ', date0         print *, 'physiq date0: ', date0
687      ENDIF test_firstcal      ENDIF test_firstcal
# Line 740  contains Line 741  contains
741    
742      ! Incrémenter le compteur de la physique      ! Incrémenter le compteur de la physique
743      itap = itap + 1      itap = itap + 1
744      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(dayvrai, 360)
745      if (julien == 0) julien = 360      if (julien == 0) julien = 360
746    
747      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
748    
749      ! Prescrire l'ozone :      ! Prescrire l'ozone :
750      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
# Line 770  contains Line 771  contains
771      frugs = MAX(frugs, 0.000015)      frugs = MAX(frugs, 0.000015)
772      zxrugs = sum(frugs * pctsrf, dim = 2)      zxrugs = sum(frugs * pctsrf, dim = 2)
773    
774      ! Calculs nécessaires au calcul de l'albedo dans l'interface      ! Calculs nécessaires au calcul de l'albedo dans l'interface avec
775        ! la surface.
776    
777      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), longi, dist)
778      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
779         CALL zenang(zlongi, time, dtphys * REAL(radpas), rmu0, fract)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
780      ELSE      ELSE
781         rmu0 = -999.999         mu0 = -999.999
782      ENDIF      ENDIF
783    
784      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
# Line 797  contains Line 799  contains
799      ! Couche limite:      ! Couche limite:
800    
801      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
802           v_seri, julien, rmu0, co2_ppm, ftsol, cdmmax, cdhmax, &           v_seri, julien, mu0, co2_ppm, ftsol, cdmmax, cdhmax, &
803           ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &           ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
804           fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &           fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &
805           fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &           fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &
# Line 931  contains Line 933  contains
933         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
934      ENDDO      ENDDO
935    
     ! Appeler la convection (au choix)  
   
     DO k = 1, llm  
        DO i = 1, klon  
           conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k) / dtphys  
           conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k) / dtphys  
        ENDDO  
     ENDDO  
   
936      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
937    
938        ! Appeler la convection (au choix)
939    
940      if (iflag_con == 2) then      if (iflag_con == 2) then
941           conv_q = d_q_dyn + d_q_vdf / dtphys
942           conv_t = d_t_dyn + d_t_vdf / dtphys
943         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
944         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
945              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
# Line 967  contains Line 964  contains
964         mfu = upwd + dnwd         mfu = upwd + dnwd
965         IF (.NOT. ok_gust) wd = 0.         IF (.NOT. ok_gust) wd = 0.
966    
967         ! Calcul des propri\'et\'es des nuages convectifs         IF (thermcep) THEN
968              zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
969         DO k = 1, llm            zqsat = zqsat / (1. - retv * zqsat)
970            DO i = 1, klon         ELSE
971               IF (thermcep) THEN            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
972                  zdelta = MAX(0., SIGN(1., rtt - t_seri(i, k)))         ENDIF
                 zqsat(i, k) = r2es * FOEEW(t_seri(i, k), zdelta) / play(i, k)  
                 zqsat(i, k) = MIN(0.5, zqsat(i, k))  
                 zqsat(i, k) = zqsat(i, k) / (1.-retv*zqsat(i, k))  
              ELSE  
                 IF (t_seri(i, k) < t_coup) THEN  
                    zqsat(i, k) = qsats(t_seri(i, k))/play(i, k)  
                 ELSE  
                    zqsat(i, k) = qsatl(t_seri(i, k))/play(i, k)  
                 ENDIF  
              ENDIF  
           ENDDO  
        ENDDO  
973    
974         ! calcul des proprietes des nuages convectifs         ! Properties of convective clouds
975         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
976         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
977              rnebcon0)              rnebcon0)
# Line 1225  contains Line 1210  contains
1210         DO i = 1, klon         DO i = 1, klon
1211            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1212            IF (thermcep) THEN            IF (thermcep) THEN
1213               zdelta = MAX(0., SIGN(1., rtt-zx_t))               zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
              zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)  
1214               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1215               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1.-retv*zx_qs)
1216               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
# Line 1245  contains Line 1229  contains
1229      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
1230      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1231         ! Get sulfate aerosol distribution :         ! Get sulfate aerosol distribution :
1232         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(dayvrai, time, firstcal, sulfate)
1233         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1234    
1235         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1236              aerindex)              aerindex)
# Line 1268  contains Line 1252  contains
1252              bl95_b1, cldtaupi, re, fl)              bl95_b1, cldtaupi, re, fl)
1253      endif      endif
1254    
     ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.  
1255      IF (MOD(itaprad, radpas) == 0) THEN      IF (MOD(itaprad, radpas) == 0) THEN
1256           ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1257         DO i = 1, klon         DO i = 1, klon
1258            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
1259                 + falbe(i, is_lic) * pctsrf(i, is_lic) &                 + falbe(i, is_lic) * pctsrf(i, is_lic) &
# Line 1281  contains Line 1265  contains
1265                 + falblw(i, is_sic) * pctsrf(i, is_sic)                 + falblw(i, is_sic) * pctsrf(i, is_sic)
1266         ENDDO         ENDDO
1267         ! Rayonnement (compatible Arpege-IFS) :         ! Rayonnement (compatible Arpege-IFS) :
1268         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, &
1269              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &
1270              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &
1271              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &
# Line 1289  contains Line 1273  contains
1273              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)
1274         itaprad = 0         itaprad = 0
1275      ENDIF      ENDIF
1276    
1277      itaprad = itaprad + 1      itaprad = itaprad + 1
1278    
1279      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
# Line 1402  contains Line 1387  contains
1387           d_qt, d_ec)           d_qt, d_ec)
1388    
1389      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1390      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, u, t, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, t, &
1391           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1392           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, albsol, rhcl, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, da, phi, mp, &
1393           cldfra, rneb, diafra, cldliq, pmflxr, pmflxs, prfl, psfl, da, phi, &           upwd, dnwd, tr_seri, zmasse)
          mp, upwd, dnwd, tr_seri, zmasse)  
1394    
1395      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &      IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1396           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, &
Legend:
Removed from v.101  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21