/[lmdze]/trunk/phylmd/phytrac.f

Diff of /trunk/phylmd/phytrac.f

Parent Directory | Revision Log | View Patch Patch

-trunk/phylmd/phytrac.f
revision 118 by guez,
Thu Dec 18 17:30:24 2014 UTC
+trunk/Sources/phylmd/phytrac.f
revision 244 by guez,
Tue Nov 14 14:56:42 2017 UTC
 Line 7 
 module phytrac_m
  contains
-   SUBROUTINE phytrac(itap, lmt_pas, julien, gmtime, firstcal, lafin, pdtphys, &
+   SUBROUTINE phytrac(julien, gmtime, firstcal, lafin, pdtphys, t_seri, paprs, &
-        t_seri, paprs, pplay, pmfu, pmfd, pde_u, pen_d, coefh, fm_therm, &
+        pplay, pmfu, pmfd, pde_u, pen_d, coefh, cdragh, fm_therm, entr_therm, &
-        entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, da, &
+        yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
-        phi, mp, upwd, dnwd, tr_seri, zmasse)
+        tr_seri, zmasse, ncid_startphy)
-     ! From phylmd/phytrac.F, version 1.15 2006/02/21 08:08:30 (SVN revision 679)
+     ! From phylmd/phytrac.F, version 1.15, 2006/02/21 08:08:30 (SVN
+     ! revision 679) and phylmd/write_histrac.h, version 1.9,
+     ! 2006/02/21 08:08:30
      ! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
-     ! Foujols, Olivia
+     ! Foujols
      ! Objet : moniteur g\'en\'eral des tendances des traceurs
-Line 23 
 contains
+Line 25 
 contains
      ! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
      ! Modifications pour les traceurs :
-     ! - uniformisation des parametrisations dans phytrac
+     ! - uniformisation des param\'etrisations dans phytrac
      ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
      use abort_gcm_m, only: abort_gcm
-     use clesphys, only: ecrit_tra
+     use clesphys2, only: conv_emanuel
-     use clesphys2, only: iflag_con
+     use cltrac_m, only: cltrac
      use cltracrn_m, only: cltracrn
+     USE conf_gcm_m, ONLY: lmt_pas
      use ctherm, only: iflag_thermals
+     use cvltr_m, only: cvltr
      use dimens_m, only: llm, nqmx
      use dimphy, only: klon
+     use histwrite_phy_m, only: histwrite_phy
      use indicesol, only: nbsrf
-     use ini_histrac_m, only: ini_histrac
+     use iniadvtrac_m, only: tname
      use initrrnpb_m, only: initrrnpb
      use minmaxqfi_m, only: minmaxqfi
+     use netcdf, only: NF90_FILL_float
+     use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
      use nflxtr_m, only: nflxtr
      use nr_util, only: assert
      use o3_chem_m, only: o3_chem
      use phyetat0_m, only: rlat
+     use phyredem0_m, only: ncid_restartphy
      use press_coefoz_m, only: press_coefoz
      use radiornpb_m, only: radiornpb
      use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
      use SUPHEC_M, only: rg
+     use time_phylmdz, only: itap
-     integer, intent(in):: itap ! number of calls to "physiq"
-     integer, intent(in):: lmt_pas ! number of time steps of "physics" per day
      integer, intent(in):: julien !jour julien, 1 <= julien <= 360
      real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
      logical, intent(in):: firstcal ! first call to "calfis"
-Line 70 
 contains
+Line 77 
 contains
      REAL pde_u(klon, llm) ! flux detraine dans le panache montant
      REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
-     REAL coefh(klon, llm) ! coeff melange couche limite
+     REAL coefh(:, 2:) ! (klon, 2:llm) coeff melange couche limite
+     real cdragh(:) ! (klon)
-     ! thermiques:
+     real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
-     real fm_therm(klon, llm+1), entr_therm(klon, llm)
+     REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau
-     ! Couche limite:
-     REAL yu1(klon) ! vents au premier niveau
-     REAL yv1(klon) ! vents au premier niveau
      ! Arguments n\'ecessaires pour les sources et puits de traceur :
-     real, intent(in):: ftsol(klon, nbsrf) ! Temperature du sol (surf)(Kelvin)
+     real, intent(in):: ftsol(:, :) ! (klon, nbsrf) surface temperature (K)
-     real pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
+     real, intent(in):: pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
      ! Lessivage pour le on-line
-     REAL frac_impa(klon, llm) ! fraction d'aerosols impactes
+     REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
-     REAL frac_nucl(klon, llm) ! fraction d'aerosols nuclees
+     REAL, intent(in):: frac_nucl(klon, llm) ! fraction d'aerosols nuclees
-     real, intent(in):: pphis(klon)
      ! Kerry Emanuel
      real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
-Line 100 
 contains
+Line 101 
 contains
      real, intent(in):: zmasse(:, :) ! (klon, llm)
      ! (column-density of mass of air in a cell, in kg m-2)
-     ! Variables local to the procedure:
+     integer, intent(in):: ncid_startphy
+     ! Local:
      integer nsplit
-Line 114 
 contains
+Line 117 
 contains
      !
      ! Pour la source de radon et son reservoir de sol
-     REAL, save:: trs(klon, nqmx - 2) ! Concentration de radon dans le sol
+     REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
      REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
      ! Masque de l'echange avec la surface
-Line 132 
 contains
+Line 135 
 contains
      SAVE scavtr
      CHARACTER itn
-     INTEGER, save:: nid_tra
-     ! nature du traceur
-     logical aerosol(nqmx - 2) ! Nature du traceur
+     logical, save:: aerosol(nqmx - 2) ! Nature du traceur
      ! ! aerosol(it) = true => aerosol
      ! ! aerosol(it) = false => gaz
-     logical clsol(nqmx - 2) ! couche limite sol calcul\'ee
-     logical radio(nqmx - 2) ! d\'ecroisssance radioactive
+     logical, save:: clsol(nqmx - 2) ! couche limite sol flux
-     save aerosol, clsol, radio
+     ! calcul\'ee, sinon prescrit
+     logical, save:: radio(nqmx - 2) ! d\'ecroisssance radioactive
      ! convection tiedtke
      INTEGER i, k, it
-Line 153 
 contains
+Line 154 
 contains
      REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
      REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
-     REAL d_tr_cv(klon, llm, nqmx - 2) ! tendance de traceurs conv pour chq traceur
+     REAL d_tr_cv(klon, llm, nqmx - 2)
+     ! tendance de traceurs conv pour chq traceur
      REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
      REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
      ! ! radioactive du rn - > pb
-     REAL d_tr_lessi_impa(klon, llm, nqmx - 2) ! la tendance du lessivage
-     ! ! par impaction
+     REAL d_tr_lessi_impa(klon, llm, nqmx - 2)
-     REAL d_tr_lessi_nucl(klon, llm, nqmx - 2) ! la tendance du lessivage
+     ! tendance du lessivage par impaction
-     ! ! par nucleation
-     REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage
+     REAL d_tr_lessi_nucl(klon, llm, nqmx - 2)
-     ! ! dans chaque couche
+     ! tendance du lessivage par nucleation
+     REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage dans chaque couche
      real ztra_th(klon, llm)
-     integer isplit
+     integer isplit, varid
      ! Controls:
-     logical:: couchelimite = .true.
      logical:: convection = .true.
-     logical:: lessivage = .true.
-     logical, save:: inirnpb
      !--------------------------------------
-Line 179 
 contains
+Line 182 
 contains
      call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
      if (firstcal) then
-        print *, 'phytrac: pdtphys = ', pdtphys
-        PRINT *, 'Frequency of tracer output: ecrit_tra = ', ecrit_tra
-        inirnpb = .true.
-        ! Initialisation des sorties :
-        call ini_histrac(nid_tra, pdtphys, nqmx - 2, lessivage)
         ! Initialisation de certaines variables pour le radon et le plomb
         ! Initialisation du traceur dans le sol (couche limite radonique)
-        trs(:, :) = 0.
+        trs(:, 2:) = 0.
-        open (unit=99, file='starttrac', status='old', err=999, &
+        call nf95_inq_varid(ncid_startphy, "trs", varid)
-             form='formatted')
+        call nf95_get_var(ncid_startphy, varid, trs(:, 1))
-        read(unit=99, fmt=*) (trs(i, 1), i=1, klon)
+        if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
-   continue
+             "some missing values in trs(:, 1)")
-        close(unit=99)
         ! Initialisation de la fraction d'aerosols lessivee
-Line 203 
 contains
+Line 198 
 contains
         ! Initialisation de la nature des traceurs
-        DO it = 1, nqmx - 2
+        aerosol = .FALSE. ! Tous les traceurs sont des gaz par defaut
-           aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
+        radio = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
-           radio(it) = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
-           clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
-        ENDDO
         if (nqmx >= 5) then
            call press_coefoz ! read input pressure levels for ozone coefficients
         end if
-     ENDIF
-     if (inirnpb) THEN
         ! Initialisation du traceur dans le sol (couche limite radonique)
         radio(1)= .true.
         radio(2)= .true.
-        clsol(1)= .true.
+        clsol(:2)= .true.
-        clsol(2)= .true.
+        clsol(3:)= .false.
         aerosol(2) = .TRUE. ! le Pb est un aerosol
         call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
-        inirnpb=.false.
      endif
      if (convection) then
         ! Calcul de l'effet de la convection
         DO it=1, nqmx - 2
-           if (iflag_con == 2) then
+           if (conv_emanuel) then
-              ! Tiedke
+              call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
+                   dnwd, d_tr_cv(:, :, it))
+           else
               CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
                    tr_seri(:, :, it), d_tr_cv(:, :, it))
-           else if (iflag_con == 3) then
-              ! Emanuel
-              call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(1, 1, it), upwd, &
-                   dnwd, d_tr_cv(1, 1, it))
            endif
            DO k = 1, llm
-Line 283 
 contains
+Line 270 
 contains
      ! Calcul de l'effet de la couche limite
-     if (couchelimite) then
+     DO k = 1, llm
-        DO k = 1, llm
+        DO i = 1, klon
-           DO i = 1, klon
+           delp(i, k) = paprs(i, k)-paprs(i, k+1)
-              delp(i, k) = paprs(i, k)-paprs(i, k+1)
-           ENDDO
         ENDDO
+     ENDDO
-        ! MAF modif pour tenir compte du cas traceur
+     ! MAF modif pour tenir compte du cas traceur
-        DO it=1, nqmx - 2
+     DO it=1, nqmx - 2
-           if (clsol(it)) then
+        if (clsol(it)) then
-              ! couche limite avec quantite dans le sol calculee
+           ! couche limite avec quantite dans le sol calculee
-              CALL cltracrn(it, pdtphys, yu1, yv1, &
+           CALL cltracrn(it, pdtphys, yu1, yv1, coefh, cdragh, t_seri, ftsol, &
-                   coefh, t_seri, ftsol, pctsrf, &
+                pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
-                   tr_seri(:, :, it), trs(1, it), &
+                masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
-                   paprs, pplay, delp, &
+                vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
-                   masktr(1, it), fshtr(1, it), hsoltr(it), &
+           DO k = 1, llm
-                   tautr(it), vdeptr(it), &
+              DO i = 1, klon
-                   rlat, &
+                 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
-                   d_tr_cl(1, 1, it), d_trs)
-              DO k = 1, llm
-                 DO i = 1, klon
-                    tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
-                 ENDDO
               ENDDO
+           ENDDO
-              ! Traceur ds sol
+           trs(:, it) = trs(:, it) + d_trs
+        else
+           ! couche limite avec flux prescrit
+           !MAF provisoire source / traceur a creer
+           DO i=1, klon
+              source(i) = 0. ! pas de source, pour l'instant
+           ENDDO
+           CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
+                paprs, pplay, delp, d_tr_cl(1, 1, it))
+           DO k = 1, llm
               DO i = 1, klon
-                 trs(i, it) = trs(i, it) + d_trs(i)
+                 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
-              END DO
-           else ! couche limite avec flux prescrit
-              !MAF provisoire source / traceur a creer
-              DO i=1, klon
-                 source(i) = 0.0 ! pas de source, pour l'instant
-              ENDDO
-              CALL cltrac(pdtphys, coefh, t_seri, &
-                   tr_seri(1, 1, it), source, &
-                   paprs, pplay, delp, &
-                   d_tr_cl(1, 1, it))
-              DO k = 1, llm
-                 DO i = 1, klon
-                    tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
-                 ENDDO
               ENDDO
-           endif
+           ENDDO
-        ENDDO
+        endif
-     endif
+     ENDDO
      ! Calcul de l'effet du puits radioactif
-Line 349 
 contains
+Line 325 
 contains
         ! Ozone as a tracer:
         if (mod(itap - 1, lmt_pas) == 0) then
            ! Once per day, update the coefficients for ozone chemistry:
-           call regr_pr_comb_coefoz(julien)
+           call regr_pr_comb_coefoz(julien, paprs, pplay)
         end if
         call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
      end if
      ! Calcul de l'effet de la precipitation
-     IF (lessivage) THEN
+     d_tr_lessi_nucl = 0.
-        d_tr_lessi_nucl = 0.
+     d_tr_lessi_impa = 0.
-        d_tr_lessi_impa = 0.
+     flestottr = 0.
-        flestottr = 0.
-        ! tendance des aerosols nuclees et impactes
-        DO it = 1, nqmx - 2
-           IF (aerosol(it)) THEN
-              DO k = 1, llm
-                 DO i = 1, klon
-                    d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
-                         (1 - frac_nucl(i, k))*tr_seri(i, k, it)
-                    d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
-                         (1 - frac_impa(i, k))*tr_seri(i, k, it)
-                 ENDDO
-              ENDDO
-           ENDIF
-        ENDDO
-        ! Mises a jour des traceurs + calcul des flux de lessivage
+     ! tendance des aerosols nuclees et impactes
-        ! Mise a jour due a l'impaction et a la nucleation
-        DO it = 1, nqmx - 2
+     DO it = 1, nqmx - 2
-           IF (aerosol(it)) THEN
+        IF (aerosol(it)) THEN
-              DO k = 1, llm
+           DO k = 1, llm
-                 DO i = 1, klon
+              DO i = 1, klon
-                    tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
+                 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
-                         * frac_nucl(i, k)
+                      (1 - frac_nucl(i, k))*tr_seri(i, k, it)
-                 ENDDO
+                 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
+                      (1 - frac_impa(i, k))*tr_seri(i, k, it)
               ENDDO
-           ENDIF
+           ENDDO
-        ENDDO
+        ENDIF
+     ENDDO
-        ! Flux lessivage total
+     ! Mises a jour des traceurs + calcul des flux de lessivage
+     ! Mise a jour due a l'impaction et a la nucleation
-        DO it = 1, nqmx - 2
+     DO it = 1, nqmx - 2
+        IF (aerosol(it)) THEN
            DO k = 1, llm
               DO i = 1, klon
-                 flestottr(i, k, it) = flestottr(i, k, it) &
+                 tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
-                      - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
+                      * frac_nucl(i, k)
-                      * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
               ENDDO
            ENDDO
+        ENDIF
+     ENDDO
+     ! Flux lessivage total
+     DO it = 1, nqmx - 2
+        DO k = 1, llm
+           DO i = 1, klon
+              flestottr(i, k, it) = flestottr(i, k, it) &
+                   - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
+                   * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
+           ENDDO
         ENDDO
-     ENDIF
+     ENDDO
      ! Ecriture des sorties
-     call write_histrac(lessivage, itap, nid_tra)
+     CALL histwrite_phy("zmasse", zmasse)
+     DO it=1, nqmx - 2
+        CALL histwrite_phy(tname(it+2), tr_seri(:, :, it))
+        CALL histwrite_phy("fl"//tname(it+2), flestottr(:, :, it))
+        CALL histwrite_phy("d_tr_th_"//tname(it+2), d_tr_th(:, :, it))
+        CALL histwrite_phy("d_tr_cv_"//tname(it+2), d_tr_cv(:, :, it))
+        CALL histwrite_phy("d_tr_cl_"//tname(it+2), d_tr_cl(:, :, it))
+     ENDDO
      if (lafin) then
-        print *, "C'est la fin de la physique."
+        call nf95_inq_varid(ncid_restartphy, "trs", varid)
-        open(unit=99, file='restarttrac', form='formatted')
+        call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
-        do i=1, klon
-           write(unit=99, fmt=*) trs(i, 1)
-        enddo
-        PRINT *, 'Ecriture du fichier restarttrac'
-        close(unit=99)
      endif
-   contains
-     subroutine write_histrac(lessivage, itap, nid_tra)
-       ! From phylmd/write_histrac.h, version 1.9 2006/02/21 08:08:30
-       use dimens_m, only: iim, jjm, llm
-       use histsync_m, only: histsync
-       use histwrite_m, only: histwrite
-       use temps, only: itau_phy
-       use iniadvtrac_m, only: tnom
-       use comgeomphy, only: airephy
-       use dimphy, only: klon
-       use grid_change, only: gr_phy_write_2d
-       use gr_phy_write_3d_m, only: gr_phy_write_3d
-       logical, intent(in):: lessivage
-       integer, intent(in):: itap ! number of calls to "physiq"
-       integer, intent(in):: nid_tra
-       ! Variables local to the procedure:
-       integer it
-       integer itau_w ! pas de temps ecriture
-       logical, parameter:: ok_sync = .true.
-       !-----------------------------------------------------
-       itau_w = itau_phy + itap
-       CALL histwrite(nid_tra, "phis", itau_w, gr_phy_write_2d(pphis))
-       CALL histwrite(nid_tra, "aire", itau_w, gr_phy_write_2d(airephy))
-       CALL histwrite(nid_tra, "zmasse", itau_w, gr_phy_write_3d(zmasse))
-       DO it=1, nqmx - 2
-          CALL histwrite(nid_tra, tnom(it+2), itau_w, &
-               gr_phy_write_3d(tr_seri(:, :, it)))
-          if (lessivage) THEN
-             CALL histwrite(nid_tra, "fl"//tnom(it+2), itau_w, &
-                  gr_phy_write_3d(flestottr(:, :, it)))
-          endif
-          CALL histwrite(nid_tra, "d_tr_th_"//tnom(it+2), itau_w, &
-               gr_phy_write_3d(d_tr_th(:, :, it)))
-          CALL histwrite(nid_tra, "d_tr_cv_"//tnom(it+2), itau_w, &
-               gr_phy_write_3d(d_tr_cv(:, :, it)))
-          CALL histwrite(nid_tra, "d_tr_cl_"//tnom(it+2), itau_w, &
-               gr_phy_write_3d(d_tr_cl(:, :, it)))
-       ENDDO
-       CALL histwrite(nid_tra, "pplay", itau_w, gr_phy_write_3d(pplay))
-       CALL histwrite(nid_tra, "T", itau_w, gr_phy_write_3d(t_seri))
-       if (ok_sync) then
-          call histsync(nid_tra)
-       endif
-     end subroutine write_histrac
    END SUBROUTINE phytrac
  end module phytrac_m

 Legend:



Removed from v.118
 


changed lines


 
Added in v.244
 Legend:



Removed from v.118
 


changed lines


 
Added in v.244
-Removed from v.118
+Added in v.244

	ViewVC Help
Powered by ViewVC 1.1.21