trunk/dyn3d/guide.f

MODULE guide_m

  ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
  ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06

  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
  REAL :: lat_min_guide, lat_max_guide

  LOGICAL :: ncep, ini_anal
  INTEGER :: online

CONTAINS

  SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)

    ! Author: F.Hourdin

    USE dimens_m, ONLY : jjm, llm
    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
    USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
    USE comvert, ONLY : ap, bp, preff, presnivs
    USE conf_gcm_m, ONLY : day_step, iperiod
    USE comgeom, ONLY : aire, rlatu, rlonv
    USE serre, ONLY : clat, clon
    USE q_sat_m, ONLY : q_sat
    USE exner_hyb_m, ONLY : exner_hyb
    USE pression_m, ONLY : pression
    USE inigrads_m, ONLY : inigrads
    use netcdf, only: nf90_nowrite, nf90_open, nf90_close

    IMPLICIT NONE

    INCLUDE 'netcdf.inc'

    !   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
    REAL :: factt, ztau(ip1jmp1)

    INTEGER, INTENT (IN) :: itau
    INTEGER :: ij, l
    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
    SAVE first
    DATA 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./

    !-----------------------------------------------------------------------

    PRINT *, 'Call sequence information: guide'

    ! calcul de l'humidite saturante

    CALL pression(ip1jmp1, ap, bp, ps, p)
    CALL massdair(p, masse)
    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

       PRINT *, 'initialisation du guide '
       CALL conf_guide
       PRINT *, 'apres conf_guide'

       file = 'guide'
       CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
            180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')

       PRINT *, &
            '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'

       IF (online==-1) RETURN
       IF (online==1) THEN

          !  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
          CALL coordij(clon, clat, ilon, ilat)
          !   aire de la maille au centre du zoom
          aire_min = aire(ilon+(ilat-1)*iip1)
          !   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)
          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
       ELSE
          alpha_t = 0.
          alpha_u = 0.
          alpha_v = 0.
          alpha_p = 0.
          !           physic=.false.
       END IF

       itau_test = 1001
       step_rea = 1
       count_no_rea = 0
       ncidpl = -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 (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
       endif

       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)
       ELSE
          status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
       END IF
       status = nf_inq_dimlen(ncidpl, rid, nlev)
       PRINT *, 'nlev', nlev
       rcod = nf90_close(ncidpl)
       !   Lecture du premier etat des reanalyses.
       CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
            masserea2, psrea2, 1, nlev)
       qrea2(:, :) = max(qrea2(:, :), 0.1)


       !   Debut de l'integration temporelle:
    END IF ! first

    ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:

    ditau = real(itau)
    dday_step = real(day_step)
    WRITE (*, *) 'ditau, dday_step'
    WRITE (*, *) ditau, dday_step
    toto = 4*ditau/dday_step
    reste = toto - aint(toto)

    IF (reste==0.) THEN
       IF (itau_test==itau) THEN
          WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
          STOP
       ELSE
          vcovrea1(:, :) = vcovrea2(:, :)
          ucovrea1(:, :) = ucovrea2(:, :)
          tetarea1(:, :) = tetarea2(:, :)
          qrea1(:, :) = qrea2(:, :)

          PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
               count_no_rea, ' non lectures'
          step_rea = step_rea + 1
          itau_test = itau
          CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
               qrea2, masserea2, psrea2, 1, nlev)
          qrea2(:, :) = max(qrea2(:, :), 0.1)
          factt = dtvr*iperiod/daysec
          ztau(:) = factt/max(alpha_t(:), 1.E-10)
          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
    !    x_gcm = a * x_gcm + (1-a) * x_reanalyses

    IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'

    ditau = real(itau)
    dday_step = real(day_step)


    tau = 4*ditau/dday_step
    tau = tau - aint(tau)

    !  ucov
    IF (guide_u) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
             ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
             IF (first .AND. ini_anal) ucov(ij, l) = a
          END DO
       END DO
    END IF

    IF (guide_t) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
             teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
             IF (first .AND. ini_anal) teta(ij, l) = a
          END DO
       END DO
    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
       CALL pression(ip1jmp1, ap, bp, ps, p)
       CALL massdair(p, masse)
    END IF


    !  q
    IF (guide_q) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
             !   hum relative en % -> hum specif
             a = qsat(ij, l)*a*0.01
             q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
             IF (first .AND. ini_anal) q(ij, l) = a
          END DO
       END DO
    END IF

    ! vcov
    IF (guide_v) THEN
       DO l = 1, llm
          DO ij = 1, ip1jm
             a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
             vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
             IF (first .AND. ini_anal) vcov(ij, l) = a
          END DO
          IF (first .AND. ini_anal) vcov(ij, l) = a
       END DO
    END IF

    first = .FALSE.

  END SUBROUTINE guide

  !=======================================================================
  SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)
    !=======================================================================

    USE dimens_m, ONLY : iim, jjm
    USE paramet_m, ONLY : iip1, jjp1
    USE comconst, ONLY : pi
    USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv
    USE serre, ONLY : clat, clon, grossismx, grossismy
    IMPLICIT NONE

    !   arguments :
    INTEGER :: type
    INTEGER :: pim, pjm
    REAL :: factt, taumin, taumax
    REAL :: dxdy_, alpha(pim, pjm)
    REAL :: dxdy_min, dxdy_max

    !  local :
    REAL :: alphamin, alphamax, gamma, xi
    SAVE gamma
    INTEGER :: i, j, ilon, ilat

    LOGICAL :: first
    SAVE first
    DATA first/ .TRUE./

    REAL :: zdx(iip1, jjp1), zdy(iip1, jjp1)

    REAL :: zlat
    REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
    COMMON /comdxdy/dxdys, dxdyu, dxdyv

    IF (first) THEN
       DO j = 2, jjm
          DO i = 2, iip1
             zdx(i, j) = 0.5*(cu_2d(i-1, j)+cu_2d(i, j))/cos(rlatu(j))
          END DO
          zdx(1, j) = zdx(iip1, j)
       END DO
       DO j = 2, jjm
          DO i = 1, iip1
             zdy(i, j) = 0.5*(cv_2d(i, j-1)+cv_2d(i, j))
          END DO
       END DO
       DO i = 1, iip1
          zdx(i, 1) = zdx(i, 2)
          zdx(i, jjp1) = zdx(i, jjm)
          zdy(i, 1) = zdy(i, 2)
          zdy(i, jjp1) = zdy(i, jjm)
       END DO
       DO j = 1, jjp1
          DO i = 1, iip1
             dxdys(i, j) = sqrt(zdx(i, j)*zdx(i, j)+zdy(i, j)*zdy(i, j))
          END DO
       END DO
       DO j = 1, jjp1
          DO i = 1, iim
             dxdyu(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
          END DO
          dxdyu(iip1, j) = dxdyu(1, j)
       END DO
       DO j = 1, jjm
          DO i = 1, iip1
             dxdyv(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
          END DO
       END DO

       CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL  ')
       CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U     ')
       CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v     ')

       !   coordonnees du centre du zoom
       CALL coordij(clon, clat, ilon, ilat)
       !   aire de la maille au centre du zoom
       dxdy_min = dxdys(ilon, ilat)
       !   dxdy maximale de la maille
       dxdy_max = 0.
       DO j = 1, jjp1
          DO i = 1, iip1
             dxdy_max = max(dxdy_max, dxdys(i, j))
          END DO
       END DO

       IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
          PRINT *, 'ATTENTION modele peu zoome'
          PRINT *, 'ATTENTION on prend une constante de guidage cste'
          gamma = 0.
       ELSE
          gamma = (dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
          PRINT *, 'gamma=', gamma
          IF (gamma<1.E-5) THEN
             PRINT *, 'gamma =', gamma, '<1e-5'
             STOP
          END IF
          PRINT *, 'gamma=', gamma
          gamma = log(0.5)/log(gamma)
       END IF
    END IF

    alphamin = factt/taumax
    alphamax = factt/taumin

    DO j = 1, pjm
       DO i = 1, pim
          IF (type==1) THEN
             dxdy_ = dxdys(i, j)
             zlat = rlatu(j)*180./pi
          ELSE IF (type==2) THEN
             dxdy_ = dxdyu(i, j)
             zlat = rlatu(j)*180./pi
          ELSE IF (type==3) THEN
             dxdy_ = dxdyv(i, j)
             zlat = rlatv(j)*180./pi
          END IF
          IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
             !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
             alpha(i, j) = alphamin
          ELSE
             xi = ((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
             xi = min(xi, 1.)
             IF (lat_min_guide<=zlat .AND. zlat<=lat_max_guide) THEN
                alpha(i, j) = xi*alphamin + (1.-xi)*alphamax
             ELSE
                alpha(i, j) = 0.
             END IF
          END IF
       END DO
    END DO


    RETURN
  END SUBROUTINE tau2alpha

END MODULE guide_m
1	guez	20	MODULE guide_m
2	guez	3
3	guez	29	! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4			! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5	guez	3
6	guez	20	REAL :: tau_min_u, tau_max_u
7			REAL :: tau_min_v, tau_max_v
8			REAL :: tau_min_t, tau_max_t
9			REAL :: tau_min_q, tau_max_q
10			REAL :: tau_min_p, tau_max_p
11			REAL :: aire_min, aire_max
12	guez	3
13
14	guez	20	LOGICAL :: guide_u, guide_v, guide_t, guide_q, guide_p
15			REAL :: lat_min_guide, lat_max_guide
16	guez	3
17	guez	20	LOGICAL :: ncep, ini_anal
18			INTEGER :: online
19	guez	3
20	guez	20	CONTAINS
21	guez	3
22	guez	29	SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
23	guez	3
24	guez	29	! Author: F.Hourdin
25	guez	3
26	guez	29	USE dimens_m, ONLY : jjm, llm
27			USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
28			USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
29			USE comvert, ONLY : ap, bp, preff, presnivs
30			USE conf_gcm_m, ONLY : day_step, iperiod
31			USE comgeom, ONLY : aire, rlatu, rlonv
32			USE serre, ONLY : clat, clon
33			USE q_sat_m, ONLY : q_sat
34			USE exner_hyb_m, ONLY : exner_hyb
35			USE pression_m, ONLY : pression
36			USE inigrads_m, ONLY : inigrads
37			use netcdf, only: nf90_nowrite, nf90_open, nf90_close
38	guez	3
39	guez	29	IMPLICIT NONE
40	guez	3
41	guez	29	INCLUDE 'netcdf.inc'
42	guez	3
43	guez	29	! variables dynamiques
44			REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
45			REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
46			REAL :: q(ip1jmp1, llm) ! temperature potentielle
47			REAL :: ps(ip1jmp1) ! pression au sol
48			REAL :: masse(ip1jmp1, llm) ! masse d'air
49	guez	3
50	guez	29	! common passe pour des sorties
51			REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
52			COMMON /comdxdy/dxdys, dxdyu, dxdyv
53	guez	3
54	guez	29	! variables dynamiques pour les reanalyses.
55			REAL :: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
56			REAL :: tetarea1(ip1jmp1, llm) ! temp pot reales
57			REAL :: qrea1(ip1jmp1, llm) ! temp pot reales
58			REAL :: psrea1(ip1jmp1) ! ps
59			REAL :: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
60			REAL :: tetarea2(ip1jmp1, llm) ! temp pot reales
61			REAL :: qrea2(ip1jmp1, llm) ! temp pot reales
62			REAL :: masserea2(ip1jmp1, llm) ! masse
63			REAL :: psrea2(ip1jmp1) ! ps
64	guez	3
65	guez	29	REAL :: alpha_q(ip1jmp1)
66			REAL :: alpha_t(ip1jmp1), alpha_p(ip1jmp1)
67			REAL :: alpha_u(ip1jmp1), alpha_v(ip1jm)
68			REAL :: dday_step, toto, reste, itau_test
69			INTEGER :: step_rea, count_no_rea
70	guez	3
71	guez	29	INTEGER :: ilon, ilat
72			REAL :: factt, ztau(ip1jmp1)
73	guez	3
74	guez	29	INTEGER, INTENT (IN) :: itau
75			INTEGER :: ij, l
76			INTEGER :: ncidpl, varidpl, nlev, status
77			INTEGER :: rcod, rid
78			REAL :: ditau, tau, a
79			SAVE nlev
80	guez	3
81	guez	29	! TEST SUR QSAT
82			REAL :: p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
83			REAL :: pkf(ip1jmp1, llm)
84			REAL :: pres(ip1jmp1, llm)
85	guez	3
86	guez	29	REAL :: qsat(ip1jmp1, llm)
87			REAL :: unskap
88			REAL :: tnat(ip1jmp1, llm)
89	guez	3
90
91	guez	29	LOGICAL :: first
92			SAVE first
93			DATA first/ .TRUE./
94	guez	3
95	guez	29	SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
96			SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
97	guez	3
98	guez	29	SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
99			SAVE step_rea, count_no_rea
100	guez	3
101	guez	29	CHARACTER (10) :: file
102			INTEGER :: igrads
103			REAL :: dtgrads
104			SAVE igrads, dtgrads
105			DATA igrads, dtgrads/2, 100./
106	guez	3
107	guez	29	!-----------------------------------------------------------------------
108	guez	3
109	guez	29	PRINT *, 'Call sequence information: guide'
110	guez	3
111	guez	29	! calcul de l'humidite saturante
112	guez	3
113	guez	29	CALL pression(ip1jmp1, ap, bp, ps, p)
114			CALL massdair(p, masse)
115			PRINT *, 'OK1'
116			CALL exner_hyb(ps, p, pks, pk, pkf)
117			PRINT *, 'OK2'
118			tnat(:, :) = pk(:, :)*teta(:, :)/cpp
119			PRINT *, 'OK3'
120			unskap = 1./kappa
121			pres(:, :) = preff(pk(:, :)/cpp)*unskap
122			PRINT *, 'OK4'
123			qsat = q_sat(tnat, pres)
124	guez	3
125	guez	29	! initialisations pour la lecture des reanalyses.
126			! alpha determine la part des injections de donnees a chaque etape
127			! alpha=1 signifie pas d'injection
128			! alpha=0 signifie injection totale
129	guez	3
130	guez	29	PRINT *, 'ONLINE=', online
131			IF (online==-1) THEN
132			RETURN
133			END IF
134	guez	3
135	guez	29	IF (first) THEN
136	guez	3
137	guez	29	PRINT *, 'initialisation du guide '
138			CALL conf_guide
139			PRINT *, 'apres conf_guide'
140	guez	3
141	guez	29	file = 'guide'
142			CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
143			180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
144	guez	3
145	guez	29	PRINT *, &
146			'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
147	guez	3
148	guez	29	IF (online==-1) RETURN
149			IF (online==1) THEN
150	guez	3
151	guez	29	! Constantes de temps de rappel en jour
152			! 0.1 c'est en gros 2h30.
153			! 1e10 est une constante infinie donc en gros pas de guidage
154	guez	3
155	guez	29	! coordonnees du centre du zoom
156			CALL coordij(clon, clat, ilon, ilat)
157			! aire de la maille au centre du zoom
158			aire_min = aire(ilon+(ilat-1)*iip1)
159			! aire maximale de la maille
160			aire_max = 0.
161			DO ij = 1, ip1jmp1
162			aire_max = max(aire_max, aire(ij))
163			END DO
164			! factt = pas de temps en fraction de jour
165			factt = dtvr*iperiod/daysec
166	guez	3
167	guez	29	CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
168			CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
169			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
170			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
171			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
172	guez	3
173	guez	29	CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
174			CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
175			CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
176	guez	3
177	guez	29	! Cas ou on force exactement par les variables analysees
178			ELSE
179			alpha_t = 0.
180			alpha_u = 0.
181			alpha_v = 0.
182			alpha_p = 0.
183			! physic=.false.
184			END IF
185	guez	3
186	guez	29	itau_test = 1001
187			step_rea = 1
188			count_no_rea = 0
189			ncidpl = -99
190	guez	3
191	guez	29	! itau_test montre si l'importation a deja ete faite au rang itau
192			! lecture d'un fichier netcdf pour determiner le nombre de niveaux
193			if (guide_u) then
194			if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
195			endif
196	guez	3
197	guez	29	if (guide_v) then
198			if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
199			endif
200	guez	3
201	guez	29	if (guide_T) then
202			if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
203			endif
204	guez	3
205	guez	29	if (guide_Q) then
206			if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
207			endif
208	guez	3
209	guez	29	IF (ncep) THEN
210			status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
211			ELSE
212			status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
213			END IF
214			status = nf_inq_dimlen(ncidpl, rid, nlev)
215			PRINT *, 'nlev', nlev
216			rcod = nf90_close(ncidpl)
217			! Lecture du premier etat des reanalyses.
218			CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
219			masserea2, psrea2, 1, nlev)
220			qrea2(:, :) = max(qrea2(:, :), 0.1)
221	guez	3
222
223	guez	29	! Debut de l'integration temporelle:
224			END IF ! first
225	guez	3
226	guez	29	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
227	guez	3
228	guez	29	ditau = real(itau)
229			dday_step = real(day_step)
230			WRITE (, ) 'ditau, dday_step'
231			WRITE (, ) ditau, dday_step
232			toto = 4*ditau/dday_step
233			reste = toto - aint(toto)
234	guez	3
235	guez	29	IF (reste==0.) THEN
236			IF (itau_test==itau) THEN
237			WRITE (, ) 'deuxieme passage de advreel a itau=', itau
238			STOP
239			ELSE
240			vcovrea1(:, :) = vcovrea2(:, :)
241			ucovrea1(:, :) = ucovrea2(:, :)
242			tetarea1(:, :) = tetarea2(:, :)
243			qrea1(:, :) = qrea2(:, :)
244	guez	3
245	guez	29	PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
246			count_no_rea, ' non lectures'
247			step_rea = step_rea + 1
248			itau_test = itau
249			CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
250			qrea2, masserea2, psrea2, 1, nlev)
251			qrea2(:, :) = max(qrea2(:, :), 0.1)
252			factt = dtvr*iperiod/daysec
253			ztau(:) = factt/max(alpha_t(:), 1.E-10)
254			CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
255			CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
256			CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
257			CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
258			CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
259			CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
260			CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
261			CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
262			CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
263			CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
264			CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
265	guez	3
266	guez	29	CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
267	guez	3
268	guez	29	END IF
269			ELSE
270			count_no_rea = count_no_rea + 1
271			END IF
272	guez	3
273	guez	29	! Guidage
274			! x_gcm = a * x_gcm + (1-a) * x_reanalyses
275	guez	3
276	guez	29	IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
277	guez	3
278	guez	29	ditau = real(itau)
279			dday_step = real(day_step)
280	guez	3
281
282	guez	29	tau = 4*ditau/dday_step
283			tau = tau - aint(tau)
284	guez	3
285	guez	29	! ucov
286			IF (guide_u) THEN
287			DO l = 1, llm
288			DO ij = 1, ip1jmp1
289			a = (1.-tau)ucovrea1(ij, l) + tauucovrea2(ij, l)
290			ucov(ij, l) = (1.-alpha_u(ij))ucov(ij, l) + alpha_u(ij)a
291			IF (first .AND. ini_anal) ucov(ij, l) = a
292			END DO
293			END DO
294			END IF
295	guez	3
296	guez	29	IF (guide_t) THEN
297			DO l = 1, llm
298			DO ij = 1, ip1jmp1
299			a = (1.-tau)tetarea1(ij, l) + tautetarea2(ij, l)
300			teta(ij, l) = (1.-alpha_t(ij))teta(ij, l) + alpha_t(ij)a
301			IF (first .AND. ini_anal) teta(ij, l) = a
302			END DO
303			END DO
304			END IF
305	guez	3
306	guez	29	! P
307			IF (guide_p) THEN
308			DO ij = 1, ip1jmp1
309			a = (1.-tau)psrea1(ij) + taupsrea2(ij)
310			ps(ij) = (1.-alpha_p(ij))ps(ij) + alpha_p(ij)a
311			IF (first .AND. ini_anal) ps(ij) = a
312			END DO
313			CALL pression(ip1jmp1, ap, bp, ps, p)
314			CALL massdair(p, masse)
315			END IF
316	guez	3
317
318	guez	29	! q
319			IF (guide_q) THEN
320			DO l = 1, llm
321			DO ij = 1, ip1jmp1
322			a = (1.-tau)qrea1(ij, l) + tauqrea2(ij, l)
323			! hum relative en % -> hum specif
324			a = qsat(ij, l)a0.01
325			q(ij, l) = (1.-alpha_q(ij))q(ij, l) + alpha_q(ij)a
326			IF (first .AND. ini_anal) q(ij, l) = a
327			END DO
328			END DO
329			END IF
330	guez	3
331	guez	29	! vcov
332			IF (guide_v) THEN
333			DO l = 1, llm
334			DO ij = 1, ip1jm
335			a = (1.-tau)vcovrea1(ij, l) + tauvcovrea2(ij, l)
336			vcov(ij, l) = (1.-alpha_v(ij))vcov(ij, l) + alpha_v(ij)a
337			IF (first .AND. ini_anal) vcov(ij, l) = a
338			END DO
339			IF (first .AND. ini_anal) vcov(ij, l) = a
340			END DO
341			END IF
342	guez	3
343	guez	29	first = .FALSE.
344	guez	3
345	guez	29	END SUBROUTINE guide
346	guez	20
347	guez	3	!=======================================================================
348	guez	20	SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)
349	guez	3	!=======================================================================
350
351	guez	20	USE dimens_m, ONLY : iim, jjm
352			USE paramet_m, ONLY : iip1, jjp1
353			USE comconst, ONLY : pi
354			USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv
355			USE serre, ONLY : clat, clon, grossismx, grossismy
356			IMPLICIT NONE
357	guez	3
358			! arguments :
359	guez	20	INTEGER :: type
360			INTEGER :: pim, pjm
361			REAL :: factt, taumin, taumax
362			REAL :: dxdy_, alpha(pim, pjm)
363			REAL :: dxdy_min, dxdy_max
364	guez	3
365			! local :
366	guez	20	REAL :: alphamin, alphamax, gamma, xi
367			SAVE gamma
368			INTEGER :: i, j, ilon, ilat
369	guez	3
370	guez	20	LOGICAL :: first
371			SAVE first
372			DATA first/ .TRUE./
373	guez	3
374	guez	20	REAL :: zdx(iip1, jjp1), zdy(iip1, jjp1)
375	guez	3
376	guez	20	REAL :: zlat
377			REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
378			COMMON /comdxdy/dxdys, dxdyu, dxdyv
379	guez	3
380	guez	20	IF (first) THEN
381			DO j = 2, jjm
382			DO i = 2, iip1
383			zdx(i, j) = 0.5*(cu_2d(i-1, j)+cu_2d(i, j))/cos(rlatu(j))
384			END DO
385			zdx(1, j) = zdx(iip1, j)
386			END DO
387			DO j = 2, jjm
388			DO i = 1, iip1
389			zdy(i, j) = 0.5*(cv_2d(i, j-1)+cv_2d(i, j))
390			END DO
391			END DO
392			DO i = 1, iip1
393			zdx(i, 1) = zdx(i, 2)
394			zdx(i, jjp1) = zdx(i, jjm)
395			zdy(i, 1) = zdy(i, 2)
396			zdy(i, jjp1) = zdy(i, jjm)
397			END DO
398			DO j = 1, jjp1
399			DO i = 1, iip1
400			dxdys(i, j) = sqrt(zdx(i, j)zdx(i, j)+zdy(i, j)zdy(i, j))
401			END DO
402			END DO
403			DO j = 1, jjp1
404			DO i = 1, iim
405			dxdyu(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
406			END DO
407			dxdyu(iip1, j) = dxdyu(1, j)
408			END DO
409			DO j = 1, jjm
410			DO i = 1, iip1
411			dxdyv(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
412			END DO
413			END DO
414	guez	3
415	guez	20	CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL ')
416			CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U ')
417			CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v ')
418	guez	3
419			! coordonnees du centre du zoom
420	guez	20	CALL coordij(clon, clat, ilon, ilat)
421	guez	3	! aire de la maille au centre du zoom
422	guez	20	dxdy_min = dxdys(ilon, ilat)
423	guez	3	! dxdy maximale de la maille
424	guez	20	dxdy_max = 0.
425			DO j = 1, jjp1
426			DO i = 1, iip1
427			dxdy_max = max(dxdy_max, dxdys(i, j))
428			END DO
429			END DO
430	guez	3
431	guez	20	IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
432			PRINT *, 'ATTENTION modele peu zoome'
433			PRINT *, 'ATTENTION on prend une constante de guidage cste'
434			gamma = 0.
435			ELSE
436			gamma = (dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
437			PRINT *, 'gamma=', gamma
438			IF (gamma<1.E-5) THEN
439			PRINT *, 'gamma =', gamma, '<1e-5'
440			STOP
441			END IF
442			PRINT *, 'gamma=', gamma
443			gamma = log(0.5)/log(gamma)
444			END IF
445			END IF
446	guez	3
447	guez	20	alphamin = factt/taumax
448			alphamax = factt/taumin
449	guez	3
450	guez	20	DO j = 1, pjm
451			DO i = 1, pim
452			IF (type==1) THEN
453			dxdy_ = dxdys(i, j)
454			zlat = rlatu(j)*180./pi
455			ELSE IF (type==2) THEN
456			dxdy_ = dxdyu(i, j)
457			zlat = rlatu(j)*180./pi
458			ELSE IF (type==3) THEN
459			dxdy_ = dxdyv(i, j)
460			zlat = rlatv(j)*180./pi
461			END IF
462			IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
463	guez	3	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
464	guez	20	alpha(i, j) = alphamin
465			ELSE
466			xi = ((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
467			xi = min(xi, 1.)
468			IF (lat_min_guide<=zlat .AND. zlat<=lat_max_guide) THEN
469			alpha(i, j) = xialphamin + (1.-xi)alphamax
470			ELSE
471			alpha(i, j) = 0.
472			END IF
473			END IF
474			END DO
475			END DO
476	guez	3
477
478	guez	20	RETURN
479			END SUBROUTINE tau2alpha
480	guez	3
481	guez	20	END MODULE guide_m