/[lmdze]/trunk/dyn3d/Guide/guide.f

Diff of /trunk/dyn3d/Guide/guide.f

Parent Directory | Revision Log | View Patch Patch

-trunk/libf/dyn3d/guide.f90
revision 39 by guez,
Tue Jan 25 15:11:05 2011 UTC
+trunk/dyn3d/guide.f
revision 108 by guez,
Tue Sep 16 14:00:41 2014 UTC
 Line 5 
 MODULE guide_m
    IMPLICIT NONE
-   REAL tau_min_u, tau_max_u
-   REAL tau_min_v, tau_max_v
-   REAL tau_min_t, tau_max_t
-   REAL tau_min_q, tau_max_q
-   REAL tau_min_p, tau_max_p
    REAL aire_min, aire_max
-   LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p
-   LOGICAL ncep, ini_anal
-   INTEGER online
  CONTAINS
-   SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
+   SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
      ! Author: F.Hourdin
-     USE comconst, ONLY : cpp, daysec, dtvr, kappa
+     USE comconst, ONLY: cpp, daysec, dtvr, kappa
-     USE comvert, ONLY : ap, bp, preff, presnivs
+     USE comgeom, ONLY: aire, rlatu, rlonv
-     USE conf_gcm_m, ONLY : day_step, iperiod
+     USE conf_gcm_m, ONLY: day_step, iperiod
-     USE comgeom, ONLY : aire, rlatu, rlonv
+     use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
-     USE dimens_m, ONLY : jjm, llm
+          ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
-     USE exner_hyb_m, ONLY : exner_hyb
+          tau_min_t, tau_max_t, tau_min_q, tau_max_q, online
-     USE inigrads_m, ONLY : inigrads
+     USE dimens_m, ONLY: iim, jjm, llm
-     use netcdf, only: nf90_nowrite, nf90_open, nf90_close
+     USE disvert_m, ONLY: ap, bp, preff, presnivs
+     USE exner_hyb_m, ONLY: exner_hyb
+     use netcdf, only: nf90_nowrite
+     use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
+          nf95_open
      use nr_util, only: pi
-     USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
+     USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
-     USE q_sat_m, ONLY : q_sat
+     USE q_sat_m, ONLY: q_sat
-     USE serre, ONLY : clat, clon
+     use read_reanalyse_m, only: read_reanalyse
-     use tau2alpha_m, only: tau2alpha
+     USE serre, ONLY: clat, clon
+     use tau2alpha_m, only: tau2alpha, dxdys
-     INCLUDE 'netcdf.inc'
+     use writefield_m, only: writefield
-     !   variables dynamiques
-     REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
-     REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
-     REAL q(ip1jmp1, llm) ! temperature potentielle
-     REAL ps(ip1jmp1) ! pression  au sol
-     REAL masse(ip1jmp1, llm) ! masse d'air
-     !   common passe pour des sorties
-     REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
-     COMMON /comdxdy/dxdys, dxdyu, dxdyv
-     !   variables dynamiques pour les reanalyses.
-     REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
-     REAL tetarea1(ip1jmp1, llm) ! temp pot  reales
-     REAL qrea1(ip1jmp1, llm) ! temp pot  reales
-     REAL psrea1(ip1jmp1) ! ps
-     REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
-     REAL tetarea2(ip1jmp1, llm) ! temp pot  reales
-     REAL qrea2(ip1jmp1, llm) ! temp pot  reales
-     REAL masserea2(ip1jmp1, llm) ! masse
-     REAL psrea2(ip1jmp1) ! ps
-     REAL alpha_q(ip1jmp1)
-     REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)
-     REAL alpha_u(ip1jmp1), alpha_v(ip1jm)
-     REAL dday_step, toto, reste, itau_test
-     INTEGER step_rea, count_no_rea
-     INTEGER ilon, ilat
+     INTEGER, INTENT(IN):: itau
-     REAL factt, ztau(ip1jmp1)
+     REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
+     REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
-     INTEGER, INTENT (IN) :: itau
+     REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
-     INTEGER ij, l
+     ! température potentielle
-     INTEGER ncidpl, varidpl, nlev, status
-     INTEGER rcod, rid
-     REAL ditau, tau, a
-     SAVE nlev
-     !  TEST SUR QSAT
-     REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
-     REAL pkf(ip1jmp1, llm)
-     REAL pres(ip1jmp1, llm)
-     REAL qsat(ip1jmp1, llm)
-     REAL unskap
-     REAL tnat(ip1jmp1, llm)
-     LOGICAL:: first = .TRUE.
-     SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
-     SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
-     SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
-     SAVE step_rea, count_no_rea
-     CHARACTER (10) file
-     INTEGER igrads
-     REAL dtgrads
-     SAVE igrads, dtgrads
-     DATA igrads, dtgrads/2, 100./
-     !-----------------------------------------------------------------------
+     REAL, intent(inout):: q(:, :, :) ! (iim + 1, jjm + 1, llm)
+     REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
-     PRINT *, 'Call sequence information: guide'
+     ! Local:
-     ! calcul de l'humidite saturante
+     ! variables dynamiques pour les réanalyses
-     forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
+     REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
-     CALL massdair(p, masse)
+     ! vents covariants reanalyses
-     PRINT *, 'OK1'
-     CALL exner_hyb(ps, p, pks, pk, pkf)
-     PRINT *, 'OK2'
-     tnat(:, :) = pk(:, :)*teta(:, :)/cpp
-     PRINT *, 'OK3'
-     unskap = 1./kappa
-     pres(:, :) = preff*(pk(:, :)/cpp)**unskap
-     PRINT *, 'OK4'
-     qsat = q_sat(tnat, pres)
-     !   initialisations pour la lecture des reanalyses.
-     !    alpha determine la part des injections de donnees a chaque etape
-     !    alpha=1 signifie pas d'injection
-     !    alpha=0 signifie injection totale
-     PRINT *, 'ONLINE=', online
-     IF (online==-1) THEN
-        RETURN
-     END IF
-     IF (first) THEN
+     REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
+     REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
-        PRINT *, 'initialisation du guide '
+     REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
-        CALL conf_guide
+     ! vents covariants reanalyses
-        PRINT *, 'apres conf_guide'
+     REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
+     REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
+     REAL, save:: masserea2(ip1jmp1, llm) ! masse
+     ! alpha determine la part des injections de donnees a chaque etape
+     ! alpha=1 signifie pas d'injection
+     ! alpha=0 signifie injection totale
+     REAL, save:: alpha_q(iim + 1, jjm + 1)
+     REAL, save:: alpha_t(iim + 1, jjm + 1)
+     REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
-        file = 'guide'
+     INTEGER, save:: step_rea, count_no_rea
-        CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
-./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
-        PRINT *, &
+     INTEGER ilon, ilat
-             '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
+     REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
+     real ztau(iim + 1, jjm + 1)
-        IF (online==-1) RETURN
+     INTEGER ij, l
-        IF (online==1) THEN
+     INTEGER ncid, dimid
+     REAL tau
+     INTEGER, SAVE:: nlev
+     ! TEST SUR QSAT
+     REAL p(iim + 1, jjm + 1, llmp1)
+     real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
+     REAL qsat(iim + 1, jjm + 1, llm)
-           !  Constantes de temps de rappel en jour
+     !-----------------------------------------------------------------------
-           !  0.1 c'est en gros 2h30.
-           !  1e10  est une constante infinie donc en gros pas de guidage
-           !   coordonnees du centre du zoom
+     !!PRINT *, 'Call sequence information: guide'
+     first_call: IF (itau == 0) THEN
+        CALL conf_guide
+        IF (online) THEN
+           ! Constantes de temps de rappel en jour
+           ! coordonnees du centre du zoom
            CALL coordij(clon, clat, ilon, ilat)
-           !   aire de la maille au centre du zoom
+           ! aire de la maille au centre du zoom
-           aire_min = aire(ilon+(ilat-1)*iip1)
+           aire_min = aire(ilon+(ilat - 1) * iip1)
-           !   aire maximale de la maille
+           ! aire maximale de la maille
            aire_max = 0.
            DO ij = 1, ip1jmp1
               aire_max = max(aire_max, aire(ij))
            END DO
-           !  factt = pas de temps en fraction de jour
-           factt = dtvr*iperiod/daysec
-           CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
+           factt = dtvr * iperiod / daysec
-           CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
-           CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
-           CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
-           CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
-           CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
-           CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')
-           CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')
-           !   Cas ou on force exactement par les variables analysees
+           CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
+           CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
+           CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
+           CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
         ELSE
+           ! Cas où on force exactement par les variables analysées
            alpha_t = 0.
            alpha_u = 0.
            alpha_v = 0.
-           alpha_p = 0.
+           alpha_q = 0.
-           !           physic=.false.
         END IF
-        itau_test = 1001
         step_rea = 1
         count_no_rea = 0
-        ncidpl = -99
+        ncid = -99
-        !    itau_test    montre si l'importation a deja ete faite au rang itau
         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
-        if (guide_u) then
+        if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncid)
-           if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
+        if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncid)
-        endif
+        if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncid)
+        if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncid)
-        if (guide_v) then
-           if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
-        endif
-        if (guide_T) then
-           if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
-        endif
-        if (guide_Q) then
-           if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
-        endif
         IF (ncep) THEN
-           status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
+           call nf95_inq_dimid(ncid, 'LEVEL', dimid)
         ELSE
-           status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
+           call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
         END IF
-        status = nf_inq_dimlen(ncidpl, rid, nlev)
+        call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
         PRINT *, 'nlev', nlev
-        rcod = nf90_close(ncidpl)
+        call nf95_close(ncid)
-        !   Lecture du premier etat des reanalyses.
+        ! Lecture du premier etat des reanalyses.
         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
-             masserea2, psrea2, 1, nlev)
+             masserea2, nlev)
-        qrea2(:, :) = max(qrea2(:, :), 0.1)
+        qrea2 = max(qrea2, 0.1)
+     END IF first_call
-        !   Debut de l'integration temporelle:
-     END IF ! first
      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
-     ditau = real(itau)
+     ! Nudging fields are given 4 times per day:
-     dday_step = real(day_step)
+     IF (mod(itau, day_step / 4) == 0) THEN
-     WRITE (*, *) 'ditau, dday_step'
+        vcovrea1 = vcovrea2
-     WRITE (*, *) ditau, dday_step
+        ucovrea1 = ucovrea2
-     toto = 4*ditau/dday_step
+        tetarea1 = tetarea2
-     reste = toto - aint(toto)
+        qrea1 = qrea2
-     IF (reste==0.) THEN
+        PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
-        IF (itau_test==itau) THEN
+             count_no_rea, ' non lectures'
-           WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
+        step_rea = step_rea + 1
-           STOP
+        CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
-        ELSE
+             masserea2, nlev)
-           vcovrea1(:, :) = vcovrea2(:, :)
+        qrea2 = max(qrea2, 0.1)
-           ucovrea1(:, :) = ucovrea2(:, :)
+        factt = dtvr * iperiod / daysec
-           tetarea1(:, :) = tetarea2(:, :)
+        ztau = factt / max(alpha_t, 1E-10)
-           qrea1(:, :) = qrea2(:, :)
+        CALL writefield("aire", aire)
+        CALL writefield("dxdys", dxdys)
-           PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
+        CALL writefield("alpha_u", alpha_u)
-                count_no_rea, ' non lectures'
+        CALL writefield("alpha_t", alpha_t)
-           step_rea = step_rea + 1
+        CALL writefield("ztau", ztau)
-           itau_test = itau
+        CALL writefield("ucov", ucov)
-           CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
+        CALL writefield("ucovrea2", ucovrea2)
-                qrea2, masserea2, psrea2, 1, nlev)
+        CALL writefield("teta", teta)
-           qrea2(:, :) = max(qrea2(:, :), 0.1)
+        CALL writefield("tetarea2", tetarea2)
-           factt = dtvr*iperiod/daysec
+        CALL writefield("qrea2", qrea2)
-           ztau(:) = factt/max(alpha_t(:), 1.E-10)
+        CALL writefield("q", q)
-           CALL wrgrads(igrads, 1, aire, 'aire      ', 'aire      ')
-           CALL wrgrads(igrads, 1, dxdys, 'dxdy      ', 'dxdy      ')
-           CALL wrgrads(igrads, 1, alpha_u, 'au        ', 'au        ')
-           CALL wrgrads(igrads, 1, alpha_t, 'at        ', 'at        ')
-           CALL wrgrads(igrads, 1, ztau, 'taut      ', 'taut      ')
-           CALL wrgrads(igrads, llm, ucov, 'u         ', 'u         ')
-           CALL wrgrads(igrads, llm, ucovrea2, 'ua        ', 'ua        ')
-           CALL wrgrads(igrads, llm, teta, 'T         ', 'T         ')
-           CALL wrgrads(igrads, llm, tetarea2, 'Ta        ', 'Ta        ')
-           CALL wrgrads(igrads, llm, qrea2, 'Qa        ', 'Qa        ')
-           CALL wrgrads(igrads, llm, q, 'Q         ', 'Q         ')
-           CALL wrgrads(igrads, llm, qsat, 'QSAT      ', 'QSAT      ')
-        END IF
      ELSE
         count_no_rea = count_no_rea + 1
      END IF
-     !   Guidage
+     ! Guidage
-     !    x_gcm = a * x_gcm + (1-a) * x_reanalyses
-     IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
-     ditau = real(itau)
+     tau = mod(real(itau) / real(day_step / 4), 1.)
-     dday_step = real(day_step)
+     ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
-     tau = 4*ditau/dday_step
-     tau = tau - aint(tau)
-     !  ucov
      IF (guide_u) THEN
-        DO l = 1, llm
+        IF (itau == 0 .AND. ini_anal) then
-           DO ij = 1, ip1jmp1
+           ucov = ucovrea1
-              a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
+        else
-              ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
+           forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
-              IF (first .AND. ini_anal) ucov(ij, l) = a
+                + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
-           END DO
+                + tau * ucovrea2(:, :, l))
-        END DO
+        end IF
      END IF
      IF (guide_t) THEN
-        DO l = 1, llm
+        IF (itau == 0 .AND. ini_anal) then
-           DO ij = 1, ip1jmp1
+           teta = tetarea1
-              a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
+        else
-              teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
+           forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
-              IF (first .AND. ini_anal) teta(ij, l) = a
+                + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
-           END DO
+                + tau * tetarea2(:, :, l))
-        END DO
+        end IF
-     END IF
-     !  P
-     IF (guide_p) THEN
-        DO ij = 1, ip1jmp1
-           a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
-           ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
-           IF (first .AND. ini_anal) ps(ij) = a
-        END DO
-        forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
-        CALL massdair(p, masse)
      END IF
-     !  q
      IF (guide_q) THEN
-        DO l = 1, llm
+        ! Calcul de l'humidité saturante :
-           DO ij = 1, ip1jmp1
+        forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
-              a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
+        CALL exner_hyb(ps, p, pks, pk)
-              !   hum relative en % -> hum specif
+        qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
-              a = qsat(ij, l)*a*0.01
-              q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
+        ! humidité relative en % -> humidité spécifique
-              IF (first .AND. ini_anal) q(ij, l) = a
+        IF (itau == 0 .AND. ini_anal) then
-           END DO
+           q = qsat * qrea1 * 0.01
-        END DO
+        else
+           forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
+                + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
+                + tau * qrea2(:, :, l)) * 0.01)
+        end IF
      END IF
-     ! vcov
      IF (guide_v) THEN
-        DO l = 1, llm
+        IF (itau == 0 .AND. ini_anal) then
-           DO ij = 1, ip1jm
+           vcov = vcovrea1
-              a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
+        else
-              vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
+           forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
-              IF (first .AND. ini_anal) vcov(ij, l) = a
+                + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
-           END DO
+                + tau * vcovrea2(:, :, l))
-           IF (first .AND. ini_anal) vcov(ij, l) = a
+        end IF
-        END DO
      END IF
-     first = .FALSE.
    END SUBROUTINE guide
  END MODULE guide_m

 Legend:



Removed from v.39
 


changed lines


 
Added in v.108
 Legend:



Removed from v.39
 


changed lines


 
Added in v.108
-Removed from v.39
+Added in v.108

	ViewVC Help
Powered by ViewVC 1.1.21