/[lmdze]/trunk/Sources/dyn3d/Guide/guide.f
ViewVC logotype

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

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

trunk/libf/dyn3d/guide.f90 revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/Sources/dyn3d/Guide/guide.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC
# Line 1  Line 1 
1  module guide_m  MODULE guide_m
2    
3    ! From dyn3d/guide.F,v 1.3 2005/05/25 13:10:09    ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4    ! and dyn3d/guide.h,v 1.1.1.1 2004/05/19 12:53:06    ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5    
6    real tau_min_u,tau_max_u    IMPLICIT NONE
   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)  
   
     use dimens_m  
     use paramet_m  
     use comconst  
     use comdissnew  
     use comvert  
     use conf_gcm_m  
     use logic  
     use comgeom  
     use serre  
     use temps  
     use tracstoke  
     use ener  
     use q_sat_m, only: q_sat  
     use exner_hyb_m, only: exner_hyb  
     use pression_m, only: pression  
     use inigrads_m, only: inigrads  
   
     IMPLICIT NONE  
   
     !      ......   Version  du 10/01/98    ..........  
   
     !             avec  coordonnees  verticales hybrides  
     !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )  
   
     !=======================================================================  
     !  
     !   Auteur:  F.Hourdin  
     !   -------  
     !  
     !   Objet:  
     !   ------  
     !  
     !   GCM LMD nouvelle grille  
     !  
     !=======================================================================  
   
     !   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv  
     !        et possibilite d'appeler une fonction f(y)  a derivee tangente  
     !        hyperbolique a la  place de la fonction a derivee sinusoidale.          
   
     !   ...  Possibilite de choisir le shema de Van-leer pour l'advection de  
     !         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .  
     !  
     !-----------------------------------------------------------------------  
     !   Declarations:  
     !   -------------  
7    
8      include "netcdf.inc"  CONTAINS
9    
10      !   variables dynamiques    SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
     REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants  
     REAL 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  
   
     !IM 180305   real aire_min,aire_max  
     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)  
     !cccccccccccccccc  
   
   
     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./  
11    
12      print *,'Call sequence information: guide'      ! Author: F.Hourdin
13    
14      !-----------------------------------------------------------------------      USE comconst, ONLY: cpp, kappa
15      ! calcul de l'humidite saturante      USE comgeom, ONLY: rlatu, rlatv
16      !-----------------------------------------------------------------------      USE conf_gcm_m, ONLY: day_step
17      CALL pression( ip1jmp1, ap, bp, ps, p )      use conf_guide_m, only: guide_u, guide_v, guide_t, guide_q, ncep, &
18      call massdair(p,masse)           ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, tau_min_t, &
19      print*,'OK1'           tau_max_t, tau_min_q, tau_max_q, online, factt
20      CALL exner_hyb(ps,p,pks,pk,pkf)      USE dimens_m, ONLY: iim, jjm, llm
21      print*,'OK2'      USE disvert_m, ONLY: ap, bp, preff, presnivs
22      tnat(:,:)=pk(:,:)*teta(:,:)/cpp      USE exner_hyb_m, ONLY: exner_hyb
23      print*,'OK3'      use init_tau2alpha_m, only: init_tau2alpha
24      unskap   = 1./ kappa      use netcdf, only: nf90_nowrite
25      pres(:,:)=preff*(pk(:,:)/cpp)**unskap      use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
26      print*,'OK4'           nf95_open
27      qsat = q_sat(tnat, pres)      use nr_util, only: pi
28        USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
29        USE q_sat_m, ONLY: q_sat
30        use read_reanalyse_m, only: read_reanalyse
31        use serre, only: grossismx, grossismy
32        use tau2alpha_m, only: tau2alpha
33        use writefield_m, only: writefield
34    
35      !-----------------------------------------------------------------------      INTEGER, INTENT(IN):: itau
36        REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
37        REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
38    
39      !-----------------------------------------------------------------------      REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
40      !   initialisations pour la lecture des reanalyses.      ! température potentielle
     !    alpha determine la part des injections de donnees a chaque etape  
     !    alpha=1 signifie pas d'injection  
     !    alpha=0 signifie injection totale  
     !-----------------------------------------------------------------------  
41    
42      print*,'ONLINE=',online      REAL, intent(inout):: q(:, :, :) ! (iim + 1, jjm + 1, llm)
43      if(online.eq.-1) then      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
        return  
     endif  
   
     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.eq.-1) return  
        if (online.eq.1) then  
   
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
           !  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  
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
           !   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))  
           enddo  
           !  factt = pas de temps en fraction de jour  
           factt=dtvr*iperiod/daysec  
   
           !     subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)  
           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   ')  
44    
45            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      ! Local:
           !   Cas ou on force exactement par les variables analysees  
        else  
           alpha_T=0.  
           alpha_u=0.  
           alpha_v=0.  
           alpha_P=0.  
           !           physic=.false.  
        endif  
   
        itau_test=1001  
        step_rea=1  
        count_no_rea=0  
        ncidpl=-99  
46    
47         !    itau_test    montre si l'importation a deja ete faite au rang itau      ! variables dynamiques pour les réanalyses
        ! lecture d'un fichier netcdf pour determiner le nombre de niveaux  
        if (guide_u) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_v) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_T) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_Q) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (ncep) then  
           status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)  
        else  
           status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)  
        endif  
        status=NF_INQ_DIMLEN(ncidpl,rid,nlev)  
        print *,'nlev', nlev  
        call ncclos(ncidpl,rcod)  
        !   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:  
        !   ----------------------------------  
48    
49      endif ! first      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
50      !      ! vents covariants reanalyses
     !-----------------------------------------------------------------------  
     !----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:  
     !-----------------------------------------------------------------------  
51    
52      ditau=real(itau)      REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
53      DDAY_step=real(day_step)      REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
54      write(*,*)'ditau,dday_step'  
55      write(*,*)ditau,dday_step      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
56      toto=4*ditau/dday_step      ! vents covariants reanalyses
57      reste=toto-aint(toto)  
58      !     write(*,*)'toto,reste',toto,reste      REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
59        REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
60      if (reste.eq.0.) then      REAL, save:: masserea2(ip1jmp1, llm) ! masse
61         if (itau_test.eq.itau) then  
62            write(*,*)'deuxieme passage de advreel a itau=',itau      ! alpha détermine la part des injections de données à chaque étape
63            stop      ! alpha=0 signifie pas d'injection
64         else      ! alpha=1 signifie injection totale
65            vcovrea1(:,:)=vcovrea2(:,:)      REAL, save:: alpha_q(iim + 1, jjm + 1)
66            ucovrea1(:,:)=ucovrea2(:,:)      REAL, save:: alpha_t(iim + 1, jjm + 1)
67            tetarea1(:,:)=tetarea2(:,:)      REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
68            qrea1(:,:)=qrea2(:,:)  
69        INTEGER, save:: step_rea, count_no_rea
70            print*,'LECTURE REANALYSES, pas ',step_rea &  
71                 ,'apres ',count_no_rea,' non lectures'      INTEGER l
72            step_rea=step_rea+1      INTEGER ncid, dimid
73            itau_test=itau      REAL tau
74            call read_reanalyse(step_rea,ps &      INTEGER, SAVE:: nlev
75                 ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)  
76            qrea2(:,:)=max(qrea2(:,:),0.1)      ! TEST SUR QSAT
77            factt=dtvr*iperiod/daysec      REAL p(iim + 1, jjm + 1, llmp1)
78            ztau(:)=factt/max(alpha_T(:),1.e-10)      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
79            call wrgrads(igrads,1,aire   ,'aire      ','aire      ' )      REAL qsat(iim + 1, jjm + 1, llm)
80            call wrgrads(igrads,1,dxdys  ,'dxdy      ','dxdy      ' )  
81            call wrgrads(igrads,1,alpha_u,'au        ','au        ' )      REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
           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      ' )  
   
        endif  
     else  
        count_no_rea=count_no_rea+1  
     endif  
82    
83      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
     !   Guidage  
     !    x_gcm = a * x_gcm + (1-a) * x_reanalyses  
     !-----------------------------------------------------------------------  
84    
85      if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'      !!PRINT *, 'Call sequence information: guide'
86    
87        first_call: IF (itau == 0) THEN
88           IF (online) THEN
89              IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
90                 ! grille regulière
91                 if (guide_u) alpha_u = factt / tau_max_u
92                 if (guide_v) alpha_v = factt / tau_max_v
93                 if (guide_t) alpha_t = factt / tau_max_t
94                 if (guide_q) alpha_q = factt / tau_max_q
95              else
96                 call init_tau2alpha(dxdys, dxdyu, dxdyv)
97    
98      ditau=real(itau)               if (guide_u) then
99      dday_step=real(day_step)                  CALL tau2alpha(dxdyu, rlatu, tau_min_u, tau_max_u, alpha_u)
100                    CALL writefield("alpha_u", alpha_u)
101                 end if
102    
103                 if (guide_v) then
104                    CALL tau2alpha(dxdyv, rlatv, tau_min_v, tau_max_v, alpha_v)
105                    CALL writefield("alpha_v", alpha_v)
106                 end if
107    
108                 if (guide_t) then
109                    CALL tau2alpha(dxdys, rlatu, tau_min_t, tau_max_t, alpha_t)
110                    CALL writefield("alpha_t", alpha_t)
111                 end if
112    
113                 if (guide_q)  then
114                    CALL tau2alpha(dxdys, rlatu, tau_min_q, tau_max_q, alpha_q)
115                    CALL writefield("alpha_q", alpha_q)
116                 end if
117              end IF
118           ELSE
119              ! Cas où on force exactement par les variables analysées
120              if (guide_u) alpha_u = 1.
121              if (guide_v) alpha_v = 1.
122              if (guide_t) alpha_t = 1.
123              if (guide_q) alpha_q = 1.
124           END IF
125    
126           step_rea = 1
127           count_no_rea = 0
128    
129      tau=4*ditau/dday_step         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux :
     tau=tau-aint(tau)  
130    
131      !  ucov         if (guide_u) then
132      if (guide_u) then            call nf95_open('u.nc',Nf90_NOWRITe,ncid)
133         do l=1,llm         else if (guide_v) then
134            do ij=1,ip1jmp1            call nf95_open('v.nc',nf90_nowrite,ncid)
135               a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)         else if (guide_T) then
136               ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a            call nf95_open('T.nc',nf90_nowrite,ncid)
              if (first.and.ini_anal) ucov(ij,l)=a  
           enddo  
        enddo  
     endif  
   
     !  teta  
     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  
           enddo  
        enddo  
     endif  
   
     !  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  
        enddo  
        CALL pression(ip1jmp1,ap,bp,ps,p)  
        CALL massdair(p,masse)  
     endif  
   
   
     !  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  
           enddo  
        enddo  
     endif  
   
     ! 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  
           enddo  
           if (first.and.ini_anal) vcov(ij,l)=a  
        enddo  
     endif  
   
     !     call dump2d(iip1,jjp1,tetarea1,'TETA REA 1     ')  
     !     call dump2d(iip1,jjp1,tetarea2,'TETA REA 2     ')  
     !     call dump2d(iip1,jjp1,teta,'TETA           ')  
   
     first=.false.  
   
     return  
   end subroutine guide  
   
   !=======================================================================  
   subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)  
     !=======================================================================  
   
     use dimens_m  
     use paramet_m  
     use comconst, only: pi  
     use comgeom  
     use serre  
     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))  
           enddo  
           zdx(1,j)=zdx(iip1,j)  
        enddo  
        do j=2,jjm  
           do i=1,iip1  
              zdy(i,j)=0.5*(cv_2d(i,j-1)+cv_2d(i,j))  
           enddo  
        enddo  
        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)  
        enddo  
        do j=1,jjp1  
           do i=1,iip1  
              dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))  
           enddo  
        enddo  
        do j=1,jjp1  
           do i=1,iim  
              dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))  
           enddo  
           dxdyu(iip1,j)=dxdyu(1,j)  
        enddo  
        do j=1,jjm  
           do i=1,iip1  
              dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))  
           enddo  
        enddo  
   
        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))  
           enddo  
        enddo  
   
        if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then  
           print*,'ATTENTION modele peu zoome'  
           print*,'ATTENTION on prend une constante de guidage cste'  
           gamma=0.  
137         else         else
138            gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)            call nf95_open('hur.nc',nf90_nowrite, ncid)
139            print*,'gamma=',gamma         end if
140            if (gamma.lt.1.e-5) then  
141               print*,'gamma =',gamma,'<1e-5'         IF (ncep) THEN
142               stop            call nf95_inq_dimid(ncid, 'LEVEL', dimid)
143            endif         ELSE
144            print*,'gamma=',gamma            call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
145            gamma=log(0.5)/log(gamma)         END IF
146         endif         call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
147      endif         PRINT *, 'nlev = ', nlev
148           call nf95_close(ncid)
149      alphamin=factt/taumax  
150      alphamax=factt/taumin         ! Lecture du premier état des réanalyses :
151           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
152      do j=1,pjm              masserea2, nlev)
153         do i=1,pim         qrea2 = max(qrea2, 0.1)
154            if (type.eq.1) then      END IF first_call
155               dxdy_=dxdys(i,j)  
156               zlat=rlatu(j)*180./pi      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
157            elseif (type.eq.2) then  
158               dxdy_=dxdyu(i,j)      ! Nudging fields are given 4 times per day:
159               zlat=rlatu(j)*180./pi      IF (mod(itau, day_step / 4) == 0) THEN
160            elseif (type.eq.3) then         vcovrea1 = vcovrea2
161               dxdy_=dxdyv(i,j)         ucovrea1 = ucovrea2
162               zlat=rlatv(j)*180./pi         tetarea1 = tetarea2
163            endif         qrea1 = qrea2
164            if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then  
165               !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin         PRINT *, 'Lecture fichiers guidage, pas ', step_rea, 'apres ', &
166               alpha(i,j)=alphamin              count_no_rea, ' non lectures'
167            else         step_rea = step_rea + 1
168               xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma         CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
169               xi=min(xi,1.)              masserea2, nlev)
170               if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then         qrea2 = max(qrea2, 0.1)
                 alpha(i,j)=xi*alphamin+(1.-xi)*alphamax  
              else  
                 alpha(i,j)=0.  
              endif  
           endif  
        enddo  
     enddo  
171    
172           if (guide_u) then
173              CALL writefield("ucov", ucov)
174              CALL writefield("ucovrea2", ucovrea2)
175           end if
176    
177           if (guide_t) then
178              CALL writefield("teta", teta)
179              CALL writefield("tetarea2", tetarea2)
180           end if
181    
182           if (guide_q) then
183              CALL writefield("qrea2", qrea2)
184              CALL writefield("q", q)
185           end if
186        ELSE
187           count_no_rea = count_no_rea + 1
188        END IF
189    
190        ! Guidage
191    
192        tau = mod(real(itau) / real(day_step / 4), 1.)
193    
194        ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
195    
196        IF (guide_u) THEN
197           IF (itau == 0 .AND. ini_anal) then
198              ucov = ucovrea1
199           else
200              forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
201                   + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
202                   + tau * ucovrea2(:, :, l))
203           end IF
204        END IF
205    
206        IF (guide_t) THEN
207           IF (itau == 0 .AND. ini_anal) then
208              teta = tetarea1
209           else
210              forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
211                   + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
212                   + tau * tetarea2(:, :, l))
213           end IF
214        END IF
215    
216        IF (guide_q) THEN
217           ! Calcul de l'humidité saturante :
218           forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
219           CALL exner_hyb(ps, p, pks, pk)
220           qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
221    
222           ! humidité relative en % -> humidité spécifique
223           IF (itau == 0 .AND. ini_anal) then
224              q = qsat * qrea1 * 0.01
225           else
226              forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
227                   + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
228                   + tau * qrea2(:, :, l)) * 0.01)
229           end IF
230        END IF
231    
232        IF (guide_v) THEN
233           IF (itau == 0 .AND. ini_anal) then
234              vcov = vcovrea1
235           else
236              forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
237                   + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
238                   + tau * vcovrea2(:, :, l))
239           end IF
240        END IF
241    
242      return    END SUBROUTINE guide
   end subroutine tau2alpha  
243    
244  end module guide_m  END MODULE guide_m
Legend:
Removed from v.3  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21