Index: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/phylmd/physiq.F	(revision 1533)
+++ LMDZ4/branches/LMDZ4_AR5/libf/phylmd/physiq.F	(revision 1534)
@@ -42,4 +42,6 @@
       use radlwsw_m, only: radlwsw
 
+!IM stations CFMIP
+      USE CFMIP_point_locations
       IMPLICIT none
 c======================================================================
@@ -676,4 +678,5 @@
 cAA
       REAL coefh(klon,klev)     ! coef d'echange pour phytrac, valable pour 2<=k<=klev
+      REAL coefm(klon,klev)     ! coef d'echange pour U, V
       REAL u1(klon)             ! vents dans la premiere couche U
       REAL v1(klon)             ! vents dans la premiere couche V
@@ -986,4 +989,5 @@
       REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
       REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
+      REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
 c#ifdef histNMC
 cym   A voir plus tard !!!!
@@ -1016,6 +1020,4 @@
       REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
       REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
-c
-cIM 280405 END
 c
       INTEGER nhori, nvert, nvert1, nvert3
@@ -1159,4 +1161,30 @@
       REAL grain(1), gtsol(1), gt2m(1), gprw(1)
 
+cIM stations CFMIP
+      INTEGER, SAVE :: nCFMIP
+c$OMP THREADPRIVATE(nCFMIP)
+      INTEGER, PARAMETER :: npCFMIP=120
+      INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
+      REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
+c$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP) 
+      INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
+      REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
+c$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
+      INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
+c$OMP THREADPRIVATE(iGCM, jGCM)
+      logical, dimension(nfiles)            :: phys_out_filestations
+      logical, parameter :: lNMC=.FALSE.
+
+cIM betaCRF
+      REAL, SAVE :: pfree, beta_pbl, beta_free
+c$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
+      REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
+c$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
+      LOGICAL, SAVE :: mskocean_beta
+c$OMP THREADPRIVATE(mskocean_beta)
+      REAL, dimension(klon, klev) :: beta       ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
+      REAL, dimension(klon, klev) :: cldtaurad  ! epaisseur optique pour radlwsw,COSP
+      REAL, dimension(klon, klev) :: cldemirad  ! emissivite pour radlwsw,COSP 
+
 cIM for NMC files
       missing_val=nf90_fill_real
@@ -1443,7 +1471,41 @@
 
 c================================================================================
-
+cIM stations CFMIP
+      nCFMIP=npCFMIP
+      OPEN(98,file='npCFMIP_param.data',status='old',
+     $          form='formatted',err=999)
+      READ(98,*,end=998) nCFMIP
+998   CONTINUE
+      CLOSE(98)
+999   CONTINUE
+      IF(nCFMIP.GT.npCFMIP) THEN
+       print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
+       CALL abort
+      else
+       print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
+      ENDIF
+c
+      ALLOCATE(tabCFMIP(nCFMIP))
+      ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
+      ALLOCATE(tabijGCM(nCFMIP))
+      ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
+      ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
+c
+c lecture des nCFMIP stations CFMIP, de leur numero 
+c et des coordonnees geographiques lonCFMIP, latCFMIP
+c
+         CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, 
+     $lonCFMIP, latCFMIP)
+c
+c identification des
+c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
+c 2) indices points tabijGCM de la grille physique 1d sur klon points
+c 3) indices iGCM, jGCM de la grille physique 2d
+c
+         CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP,
+     $tabijGCM, lonGCM, latGCM, iGCM, jGCM)
+c
          ENDIF !debut
-
+ 
            DO i=1,klon
              rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
@@ -1483,6 +1545,4 @@
      .                   lmt_pas
 c
-cIM 030306 END
-
       capemaxcels = 't_max(X)'
       t2mincels = 't_min(X)'
@@ -1501,9 +1561,13 @@
 
 c$OMP MASTER
-       call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta, 
-     &                        ctetaSTD,dtime,ok_veget,
-     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
-     &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, 
-     &                        read_climoz, new_aod, aerosol_couple)
+       call phys_output_open(rlon,rlat,nCFMIP,tabijGCM,
+     &                       iGCM,jGCM,lonGCM,latGCM,
+     &                       jjmp1,nlevSTD,clevSTD,
+     &                       nbteta, ctetaSTD, dtime,ok_veget,
+     &                       type_ocean,iflag_pbl,ok_mensuel,ok_journe,
+     &                       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, 
+     &                       read_climoz, phys_out_filestations,
+     &                       new_aod, aerosol_couple
+     &                        )
 c$OMP END MASTER
 c$OMP BARRIER
@@ -1525,5 +1589,4 @@
 #endif
 
-cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf 
          ecrit_hf2mth = ecrit_mth/ecrit_hf
 
@@ -1538,5 +1601,4 @@
          ecrit_reg = ecrit_reg * un_jour
          ecrit_tra = ecrit_tra * un_jour
-         ecrit_ISCCP = ecrit_ISCCP * un_jour
          ecrit_LES = ecrit_LES * un_jour
 c
@@ -1544,6 +1606,4 @@
      .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
      .   ecrit_hf2mth
-cIM 030306 END
-
 
 cXXXPB Positionner date0 pour initialisation de ORCHIDEE
@@ -1603,4 +1663,36 @@
       END IF
 C$omp end single
+c
+cIM betaCRF
+      pfree=70000. !Pa
+      beta_pbl=1.
+      beta_free=1.
+      lon1_beta=-180.
+      lon2_beta=+180.
+      lat1_beta=90.
+      lat2_beta=-90.
+      mskocean_beta=.FALSE.
+
+      OPEN(99,file='beta_crf.data',status='old',
+     $          form='formatted',err=9999)
+      READ(99,*,end=9998) pfree
+      READ(99,*,end=9998) beta_pbl
+      READ(99,*,end=9998) beta_free
+      READ(99,*,end=9998) lon1_beta
+      READ(99,*,end=9998) lon2_beta
+      READ(99,*,end=9998) lat1_beta
+      READ(99,*,end=9998) lat2_beta
+      READ(99,*,end=9998) mskocean_beta
+9998  Continue
+      CLOSE(99)
+9999  Continue
+      WRITE(*,*)'pfree=',pfree
+      WRITE(*,*)'beta_pbl=',beta_pbl
+      WRITE(*,*)'beta_free=',beta_free
+      WRITE(*,*)'lon1_beta=',lon1_beta
+      WRITE(*,*)'lon2_beta=',lon2_beta
+      WRITE(*,*)'lat1_beta=',lat1_beta
+      WRITE(*,*)'lat2_beta=',lat2_beta
+      WRITE(*,*)'mskocean_beta=',mskocean_beta
       ENDIF
 !
@@ -1906,5 +1998,5 @@
      s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
      s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
-     s     coefh,     slab_wfbils,                
+     s     coefh,     coefm,     slab_wfbils,                
      d     qsol,      zq2m,      s_pblh,  s_lcl,
      d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
@@ -2933,5 +3025,5 @@
          mass_solu_aero_pi(:,:) = ccm(:,:,2) 
       END IF
-
+c
       if (ok_newmicro) then
       CALL newmicro (paprs, pplay,ok_newmicro,
@@ -2954,4 +3046,55 @@
       endif
 c
+cIM betaCRF
+c
+      cldtaurad = cldtau
+      cldemirad = cldemi
+c
+      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
+     $lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
+c
+c global
+c
+       DO k=1, klev
+       DO i=1, klon
+        if (pplay(i,k).GE.pfree) THEN
+         beta(i,k) = beta_pbl
+        else
+         beta(i,k) = beta_free
+        endif
+        if (mskocean_beta) THEN
+         beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
+        endif
+        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
+        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
+       ENDDO
+       ENDDO 
+c
+      else
+c
+c regional
+c
+       DO k=1, klev
+       DO i=1,klon
+c
+        if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND.
+     $      rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
+         if (pplay(i,k).GE.pfree) THEN
+          beta(i,k) = beta_pbl
+         else
+          beta(i,k) = beta_free
+         endif
+         if (mskocean_beta) THEN
+          beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
+         endif
+        cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
+        cldemirad(i,k) = cldemi(i,k) * beta(i,k)
+        endif
+c
+       ENDDO
+       ENDDO
+c
+      endif
+c
 c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
 c
@@ -2982,5 +3125,5 @@
      e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
      e        wo(:, :, 1),
-     e        cldfra, cldemi, cldtau,
+     e        cldfra, cldemirad, cldtaurad,
      s        heat,heat0,cool,cool0,radsol,albpla,
      s        topsw,toplw,solsw,sollw,
@@ -3000,10 +3143,18 @@
 #endif
       ELSE
-
+c
+cIM calcul radiatif pour le cas actuel
+c
+       RCO2 = RCO2_act
+       RCH4 = RCH4_act
+       RN2O = RN2O_act
+       RCFC11 = RCFC11_act
+       RCFC12 = RCFC12_act
+c
          CALL radlwsw
      e        (dist, rmu0, fract, 
      e        paprs, pplay,zxtsol,albsol1, albsol2, 
      e        t_seri,q_seri,wo,
-     e        cldfra, cldemi, cldtau,
+     e        cldfra, cldemirad, cldtaurad,
      e        ok_ade, ok_aie,
      e        tau_aero, piz_aero, cg_aero,
@@ -3023,5 +3174,42 @@
      o        topswcf_aero, solswcf_aero)
          
-
+c
+cIM 2eme calcul radiatif pour le cas perturbe ou au moins un
+cIM des taux doit etre different du taux actuel
+cIM Par defaut on a les taux perturbes egaux aux taux actuels
+c
+       if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR.
+     $RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR.
+     $RCFC12_per.NE.RCFC12_act) THEN
+c
+       RCO2 = RCO2_per
+       RCH4 = RCH4_per
+       RN2O = RN2O_per
+       RCFC11 = RCFC11_per
+       RCFC12 = RCFC12_per
+c
+         CALL radlwsw
+     e        (dist, rmu0, fract, 
+     e        paprs, pplay,zxtsol,albsol1, albsol2, 
+     e        t_seri,q_seri,wo,
+     e        cldfra, cldemi, cldtau,
+     e        ok_ade, ok_aie,
+     e        tau_aero, piz_aero, cg_aero,
+     e        cldtaupi,new_aod,
+     e        zqsat, flwc, fiwc,
+     s        heatp,heat0p,coolp,cool0p,radsolp,albplap,
+     s        topswp,toplwp,solswp,sollwp,
+     s        sollwdownp,
+     s        topsw0p,toplw0p,solsw0p,sollw0p,
+     s        lwdn0p, lwdnp, lwup0p, lwupp, 
+     s        swdn0p, swdnp, swup0p, swupp,
+     s        topswad_aerop, solswad_aerop,
+     s        topswai_aerop, solswai_aerop,
+     o        topswad0_aerop, solswad0_aerop,
+     o        topsw_aerop, topsw0_aerop,
+     o        solsw_aerop, solsw0_aerop,
+     o        topswcf_aerop, solswcf_aerop)
+       endif
+c
       ENDIF ! aerosol_couple
       itaprad = 0
@@ -3184,5 +3372,5 @@
 c
 c  ajout des tendances
-        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
+        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
 
       ENDIF
@@ -3260,5 +3448,5 @@
      $                   prfl(:,1:klev),psfl(:,1:klev),
      $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
-     $                   mr_ozone,cldtau, cldemi)
+     $                   mr_ozone,cldtaurad, cldemirad)
 
 !     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
@@ -3416,4 +3604,13 @@
       vwriteSTD(:,:,5)=vlevSTD(:,:)
       wwriteSTD(:,:,5)=wlevSTD(:,:)
+c
+cIM initialisation 6eme fichier de sortie 
+      twriteSTD(:,:,6)=tlevSTD(:,:)
+      qwriteSTD(:,:,6)=qlevSTD(:,:)
+      rhwriteSTD(:,:,6)=rhlevSTD(:,:)
+      phiwriteSTD(:,:,6)=philevSTD(:,:)
+      uwriteSTD(:,:,6)=ulevSTD(:,:)
+      vwriteSTD(:,:,6)=vlevSTD(:,:)
+      wwriteSTD(:,:,6)=wlevSTD(:,:)
 cIM for NMC files
       DO n=1, nlevSTD3
