Ignore:
Timestamp:
Jul 5, 2013, 10:38:13 AM (11 years ago)
Author:
Laurent Fairhead
Message:

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL

AMMA and Sandu cases included in 1D configuration
MPL

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phy1d/1DUTILS.h_no_writelim

    r1763 r1780  
    2727#include "compar1d.h"
    2828#include "flux_arp.h"
     29#include "tsoilnudge.h"
    2930#include "fcg_gcssold.h"
    3031#include "fcg_racmo.h"
    … …  
    100101!             initial profiles from RICO idealized
    101102!             LS convergence imposed from  RICO (cst)
     103!         = 6 ==> forcing_amma = .true.
    102104!         = 40 ==> forcing_GCSSold = .true.
    103105!             initial profile from GCSS file
    104106!             LS convergence imposed from GCSS file
     107!         = 59 ==> forcing_sandu = .true.
     108!             initial profiles from sanduref file: see prof.inp.001
     109!             SST varying with time and divergence constante: see ifa_sanduref.txt file
     110!             Radiation has to be computed interactively
     111!         = 60 ==> forcing_astex = .true.
     112!             initial profiles from file: see prof.inp.001
     113!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
     114!             Radiation has to be computed interactively
     115!         = 61 ==> forcing_armcu = .true.
     116!             initial profiles from file: see prof.inp.001
     117!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
     118!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
     119!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
     120!             Radiation to be switched off
    105121!
    106122       forcing_type = 0
    … …  
    234250       zpicinp = 300.
    235251       CALL getin('zpicinp',zpicinp)
     252!Config key = nudge_tsoil
     253!Config  Desc = activation of soil temperature nudging
     254!Config  Def  = .FALSE.
     255!Config  Help = ...
     256
     257       nudge_tsoil=.FALSE.
     258       CALL getin('nudge_tsoil',nudge_tsoil)
     259
     260!Config key = isoil_nudge
     261!Config  Desc = level number where soil temperature is nudged
     262!Config  Def  = 3
     263!Config  Help = ...
     264
     265       isoil_nudge=3
     266       CALL getin('isoil_nudge',isoil_nudge)
     267
     268!Config key = Tsoil_nudge
     269!Config  Desc = target temperature for tsoil(isoil_nudge)
     270!Config  Def  = 300.
     271!Config  Help = ...
     272
     273       Tsoil_nudge=300.
     274       CALL getin('Tsoil_nudge',Tsoil_nudge)
     275
     276!Config key = tau_soil_nudge
     277!Config  Desc = nudging relaxation time for tsoil
     278!Config  Def  = 3600.
     279!Config  Help = ...
     280
     281       tau_soil_nudge=3600.
     282       CALL getin('tau_soil_nudge',tau_soil_nudge)
     283
     284
    236285
    237286
    … …  
    257306      write(lunout,*)' qsolinp = ', qsolinp
    258307      write(lunout,*)' zpicinp = ', zpicinp
     308      write(lunout,*)' nudge_tsoil = ', nudge_tsoil
     309      write(lunout,*)' isoil_nudge = ', isoil_nudge
     310      write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
     311      write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
    259312      IF (forcing_type .eq.40) THEN
    260313        write(lunout,*) '--- Forcing type GCSS Old --- with:'
    … …  
    9651018      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    9661019 
     1020!    Ancienne version disvert dont on a modifie nom pour utiliser
     1021!    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
     1022!    (MPL 18092012)
     1023!
    9671024!    Auteur :  P. Le Van .
    9681025!
    … …  
    14091466          end
    14101467
     1468c-------------------------------------------------------------------------
     1469      SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
     1470      implicit none
     1471
     1472c-------------------------------------------------------------------------
     1473c Read I.SANDU case forcing data
     1474c-------------------------------------------------------------------------
     1475
     1476      integer nlev_sandu,nt_sandu
     1477      real ts_sandu(nt_sandu)
     1478      character*80 fich_sandu
     1479
     1480      integer no,l,k,ip
     1481      real riy,rim,rid,rih,bid
     1482
     1483      integer iy,im,id,ih
     1484
     1485       real plev_min
     1486
     1487       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
     1488
     1489      open(21,file=trim(fich_sandu),form='formatted')
     1490      read(21,'(a)')
     1491      do ip = 1, nt_sandu
     1492      read(21,'(a)')
     1493      read(21,'(a)')
     1494      read(21,223) iy, im, id, ih, ts_sandu(ip)
     1495      print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
     1496      enddo
     1497      close(21)
     1498
     1499  223 format(4i3,f8.2)
     1500  226 format(f7.1,1x,10f8.2)
     1501  227 format(f7.1,1x,1p,4e11.3)
     1502  230 format(6f9.3,4e11.3)
     1503
     1504          return
     1505          end
     1506
     1507!=====================================================================
     1508c-------------------------------------------------------------------------
     1509      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,
     1510     : ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
     1511      implicit none
     1512
     1513c-------------------------------------------------------------------------
     1514c Read Astex case forcing data
     1515c-------------------------------------------------------------------------
     1516
     1517      integer nlev_astex,nt_astex
     1518      real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
     1519      real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
     1520      character*80 fich_astex
     1521
     1522      integer no,l,k,ip
     1523      real riy,rim,rid,rih,bid
     1524
     1525      integer iy,im,id,ih
     1526
     1527       real plev_min
     1528
     1529       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
     1530
     1531      open(21,file=trim(fich_astex),form='formatted')
     1532      read(21,'(a)')
     1533      read(21,'(a)')
     1534      do ip = 1, nt_astex
     1535      read(21,'(a)')
     1536      read(21,'(a)')
     1537      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),
     1538     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
     1539      ts_astex(ip)=ts_astex(ip)+273.15
     1540      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),
     1541     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
     1542      enddo
     1543      close(21)
     1544
     1545  223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
     1546  226 format(f7.1,1x,10f8.2)
     1547  227 format(f7.1,1x,1p,4e11.3)
     1548  230 format(6f9.3,4e11.3)
     1549
     1550          return
     1551          end
    14111552!=====================================================================
    14121553      subroutine read_twpice(fich_twpice,nlevel,ntime
    … …  
    18912032       return
    18922033       end
     2034!=====================================================================
     2035
     2036       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof
     2037     :         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof
     2038     :         ,omega_prof,o3mmr_prof
     2039     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
     2040     :         ,omega_mod,o3mmr_mod,mxcalc)
     2041
     2042       implicit none
     2043
     2044#include "dimensions.h"
     2045
     2046c-------------------------------------------------------------------------
     2047c Vertical interpolation of SANDUREF forcing data onto model levels
     2048c-------------------------------------------------------------------------
     2049
     2050       integer nlevmax
     2051       parameter (nlevmax=41)
     2052       integer nlev_sandu,mxcalc
     2053!       real play(llm), plev_prof(nlevmax)
     2054!       real t_prof(nlevmax),q_prof(nlevmax)
     2055!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
     2056!       real ht_prof(nlevmax),vt_prof(nlevmax)
     2057!       real hq_prof(nlevmax),vq_prof(nlevmax)
     2058
     2059       real play(llm), plev_prof(nlev_sandu)
     2060       real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
     2061       real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
     2062       real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
     2063
     2064       real t_mod(llm),thl_mod(llm),q_mod(llm)
     2065       real u_mod(llm),v_mod(llm), w_mod(llm)
     2066       real omega_mod(llm),o3mmr_mod(llm)
     2067
     2068       integer l,k,k1,k2,kp
     2069       real aa,frac,frac1,frac2,fact
     2070
     2071       do l = 1, llm
     2072
     2073        if (play(l).ge.plev_prof(nlev_sandu)) then
     2074
     2075        mxcalc=l
     2076         k1=0
     2077         k2=0
     2078
     2079         if (play(l).le.plev_prof(1)) then
     2080
     2081         do k = 1, nlev_sandu-1
     2082          if (play(l).le.plev_prof(k)
     2083     :       .and. play(l).gt.plev_prof(k+1)) then
     2084            k1=k
     2085            k2=k+1
     2086          endif
     2087         enddo
     2088
     2089         if (k1.eq.0 .or. k2.eq.0) then
     2090          write(*,*) 'PB! k1, k2 = ',k1,k2
     2091          write(*,*) 'l,play(l) = ',l,play(l)/100
     2092         do k = 1, nlev_sandu-1
     2093          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
     2094         enddo
     2095         endif
     2096
     2097         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
     2098         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
     2099         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
     2100         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))
     2101         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
     2102         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
     2103         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
     2104         omega_mod(l)=omega_prof(k2)-
     2105     :      frac*(omega_prof(k2)-omega_prof(k1))
     2106         o3mmr_mod(l)=o3mmr_prof(k2)-
     2107     :      frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
     2108
     2109         else !play>plev_prof(1)
     2110
     2111         k1=1
     2112         k2=2
     2113         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
     2114         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
     2115         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
     2116         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
     2117         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
     2118         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
     2119         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
     2120         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
     2121         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
     2122         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
     2123
     2124         endif ! play.le.plev_prof(1)
     2125
     2126        else ! above max altitude of forcing file
     2127
     2128cjyg
     2129         fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
     2130         fact = max(fact,0.)                                           !jyg
     2131         fact = exp(-fact)                                             !jyg
     2132         t_mod(l)= t_prof(nlev_sandu)                                   !jyg
     2133         thl_mod(l)= thl_prof(nlev_sandu)                                   !jyg
     2134         q_mod(l)= q_prof(nlev_sandu)*fact                              !jyg
     2135         u_mod(l)= u_prof(nlev_sandu)*fact                              !jyg
     2136         v_mod(l)= v_prof(nlev_sandu)*fact                              !jyg
     2137         w_mod(l)= w_prof(nlev_sandu)*fact                              !jyg
     2138         omega_mod(l)= omega_prof(nlev_sandu)*fact                      !jyg
     2139         o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact                      !jyg
     2140
     2141        endif ! play
     2142
     2143       enddo ! l
     2144
     2145       do l = 1,llm
     2146!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
     2147!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
     2148       enddo
     2149
     2150          return
     2151          end
     2152!=====================================================================
     2153       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof
     2154     :         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof
     2155     :         ,w_prof,tke_prof,o3mmr_prof
     2156     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
     2157     :         ,tke_mod,o3mmr_mod,mxcalc)
     2158
     2159       implicit none
     2160
     2161#include "dimensions.h"
     2162
     2163c-------------------------------------------------------------------------
     2164c Vertical interpolation of Astex forcing data onto model levels
     2165c-------------------------------------------------------------------------
     2166
     2167       integer nlevmax
     2168       parameter (nlevmax=41)
     2169       integer nlev_astex,mxcalc
     2170!       real play(llm), plev_prof(nlevmax)
     2171!       real t_prof(nlevmax),qv_prof(nlevmax)
     2172!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
     2173!       real ht_prof(nlevmax),vt_prof(nlevmax)
     2174!       real hq_prof(nlevmax),vq_prof(nlevmax)
     2175
     2176       real play(llm), plev_prof(nlev_astex)
     2177       real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
     2178       real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
     2179       real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
     2180       real qt_prof(nlev_astex),tke_prof(nlev_astex)
     2181
     2182       real t_mod(llm),thl_mod(llm),qv_mod(llm)
     2183       real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
     2184       real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)
     2185
     2186       integer l,k,k1,k2,kp
     2187       real aa,frac,frac1,frac2,fact
     2188
     2189       do l = 1, llm
     2190
     2191        if (play(l).ge.plev_prof(nlev_astex)) then
     2192
     2193        mxcalc=l
     2194         k1=0
     2195         k2=0
     2196
     2197         if (play(l).le.plev_prof(1)) then
     2198
     2199         do k = 1, nlev_astex-1
     2200          if (play(l).le.plev_prof(k)
     2201     :       .and. play(l).gt.plev_prof(k+1)) then
     2202            k1=k
     2203            k2=k+1
     2204          endif
     2205         enddo
     2206
     2207         if (k1.eq.0 .or. k2.eq.0) then
     2208          write(*,*) 'PB! k1, k2 = ',k1,k2
     2209          write(*,*) 'l,play(l) = ',l,play(l)/100
     2210         do k = 1, nlev_astex-1
     2211          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
     2212         enddo
     2213         endif
     2214
     2215         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
     2216         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
     2217         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
     2218         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
     2219         ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
     2220         qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))
     2221         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
     2222         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
     2223         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
     2224         tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
     2225         o3mmr_mod(l)=o3mmr_prof(k2)-
     2226     :      frac*(o3mmr_prof(k2)-o3mmr_prof(k1))
     2227
     2228         else !play>plev_prof(1)
     2229
     2230         k1=1
     2231         k2=2
     2232         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
     2233         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
     2234         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
     2235         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
     2236         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
     2237         ql_mod(l)= frac1*ql_prof(k1) - frac2*ql_prof(k2)
     2238         qt_mod(l)= frac1*qt_prof(k1) - frac2*qt_prof(k2)
     2239         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
     2240         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
     2241         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
     2242         tke_mod(l)= frac1*tke_prof(k1) - frac2*tke_prof(k2)
     2243         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
     2244
     2245         endif ! play.le.plev_prof(1)
     2246
     2247        else ! above max altitude of forcing file
     2248
     2249cjyg
     2250         fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
     2251         fact = max(fact,0.)                                           !jyg
     2252         fact = exp(-fact)                                             !jyg
     2253         t_mod(l)= t_prof(nlev_astex)                                   !jyg
     2254         thl_mod(l)= thl_prof(nlev_astex)                                   !jyg
     2255         qv_mod(l)= qv_prof(nlev_astex)*fact                              !jyg
     2256         ql_mod(l)= ql_prof(nlev_astex)*fact                              !jyg
     2257         qt_mod(l)= qt_prof(nlev_astex)*fact                              !jyg
     2258         u_mod(l)= u_prof(nlev_astex)*fact                              !jyg
     2259         v_mod(l)= v_prof(nlev_astex)*fact                              !jyg
     2260         w_mod(l)= w_prof(nlev_astex)*fact                              !jyg
     2261         tke_mod(l)= tke_prof(nlev_astex)*fact                              !jyg
     2262         o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact                      !jyg
     2263
     2264        endif ! play
     2265
     2266       enddo ! l
     2267
     2268       do l = 1,llm
     2269!      print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
     2270!    $        l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
     2271       enddo
     2272
     2273          return
     2274          end
    18932275
    18942276!======================================================================
    … …  
    20552437          end
    20562438
     2439!======================================================================
     2440        SUBROUTINE interp_sandu_time(day,day1,annee_ref
     2441     i             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu
     2442     i             ,nlev_sandu,ts_sandu,ts_prof)
     2443        implicit none
     2444
     2445!---------------------------------------------------------------------------------------
     2446! Time interpolation of a 2D field to the timestep corresponding to day
     2447!
     2448! day: current julian day (e.g. 717538.2)
     2449! day1: first day of the simulation
     2450! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
     2451! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
     2452!---------------------------------------------------------------------------------------
     2453! inputs:
     2454        integer annee_ref
     2455        integer nt_sandu,nlev_sandu
     2456        integer year_ini_sandu
     2457        real day, day1,day_ini_sandu,dt_sandu
     2458        real ts_sandu(nt_sandu)
     2459! outputs:
     2460        real ts_prof
     2461! local:
     2462        integer it_sandu1, it_sandu2,k
     2463        real timeit,time_sandu1,time_sandu2,frac
     2464! Check that initial day of the simulation consistent with SANDU period:
     2465       if (annee_ref.ne.2006 ) then
     2466        print*,'Pour SANDUREF, annee_ref doit etre 2006 '
     2467        print*,'Changer annee_ref dans run.def'
     2468        stop
     2469       endif
     2470!      if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
     2471!       print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
     2472!       print*,'Changer dayref dans run.def'
     2473!       stop
     2474!      endif
     2475
     2476! Determine timestep relative to the 1st day of TOGA-COARE:
     2477!       timeit=(day-day1)*86400.
     2478!       if (annee_ref.eq.1992) then
     2479!        timeit=(day-day_ini_sandu)*86400.
     2480!       else
     2481!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
     2482!       endif
     2483      timeit=(day-day_ini_sandu)*86400
     2484
     2485! Determine the closest observation times:
     2486       it_sandu1=INT(timeit/dt_sandu)+1
     2487       it_sandu2=it_sandu1 + 1
     2488       time_sandu1=(it_sandu1-1)*dt_sandu
     2489       time_sandu2=(it_sandu2-1)*dt_sandu
     2490       print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
     2491       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',
     2492     .          it_sandu1,it_sandu2,time_sandu1,time_sandu2
     2493
     2494       if (it_sandu1 .ge. nt_sandu) then
     2495        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '
     2496     :        ,day,it_sandu1,it_sandu2,timeit/86400.
     2497        stop
     2498       endif
     2499
     2500! time interpolation:
     2501       frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
     2502       frac=max(frac,0.0)
     2503
     2504       ts_prof = ts_sandu(it_sandu2)
     2505     :          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
     2506
     2507         print*,
     2508     :'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',
     2509     :day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,
     2510     :it_sandu2,ts_prof
     2511
     2512        return
     2513        END
    20572514!=====================================================================
    20582515c-------------------------------------------------------------------------
    … …  
    22162673          end
    22172674 
     2675!======================================================================
     2676        SUBROUTINE interp_astex_time(day,day1,annee_ref
     2677     i             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex
     2678     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
     2679     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
     2680     i             ,ufa_prof,vfa_prof)
     2681        implicit none
     2682
     2683!---------------------------------------------------------------------------------------
     2684! Time interpolation of a 2D field to the timestep corresponding to day
     2685!
     2686! day: current julian day (e.g. 717538.2)
     2687! day1: first day of the simulation
     2688! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)
     2689! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
     2690!---------------------------------------------------------------------------------------
     2691
     2692! inputs:
     2693        integer annee_ref
     2694        integer nt_astex,nlev_astex
     2695        integer year_ini_astex
     2696        real day, day1,day_ini_astex,dt_astex
     2697        real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
     2698        real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
     2699! outputs:
     2700        real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
     2701! local:
     2702        integer it_astex1, it_astex2,k
     2703        real timeit,time_astex1,time_astex2,frac
     2704
     2705! Check that initial day of the simulation consistent with ASTEX period:
     2706       if (annee_ref.ne.1992 ) then
     2707        print*,'Pour Astex, annee_ref doit etre 1992 '
     2708        print*,'Changer annee_ref dans run.def'
     2709        stop
     2710       endif
     2711       if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
     2712        print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
     2713        print*,'Changer dayref dans run.def'
     2714        stop
     2715       endif
     2716
     2717! Determine timestep relative to the 1st day of TOGA-COARE:
     2718!       timeit=(day-day1)*86400.
     2719!       if (annee_ref.eq.1992) then
     2720!        timeit=(day-day_ini_astex)*86400.
     2721!       else
     2722!        timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
     2723!       endif
     2724      timeit=(day-day_ini_astex)*86400
     2725
     2726! Determine the closest observation times:
     2727       it_astex1=INT(timeit/dt_astex)+1
     2728       it_astex2=it_astex1 + 1
     2729       time_astex1=(it_astex1-1)*dt_astex
     2730       time_astex2=(it_astex2-1)*dt_astex
     2731       print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
     2732       print *,'it_astex1,it_astex2,time_astex1,time_astex2',
     2733     .          it_astex1,it_astex2,time_astex1,time_astex2
     2734
     2735       if (it_astex1 .ge. nt_astex) then
     2736        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '
     2737     :        ,day,it_astex1,it_astex2,timeit/86400.
     2738        stop
     2739       endif
     2740
     2741! time interpolation:
     2742       frac=(time_astex2-timeit)/(time_astex2-time_astex1)
     2743       frac=max(frac,0.0)
     2744
     2745       div_prof = div_astex(it_astex2)
     2746     :          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
     2747       ts_prof = ts_astex(it_astex2)
     2748     :          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
     2749       ug_prof = ug_astex(it_astex2)
     2750     :          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
     2751       vg_prof = vg_astex(it_astex2)
     2752     :          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
     2753       ufa_prof = ufa_astex(it_astex2)
     2754     :          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
     2755       vfa_prof = vfa_astex(it_astex2)
     2756     :          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
     2757
     2758         print*,
     2759     :'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',
     2760     :day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,
     2761     :it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
     2762
     2763        return
     2764        END
     2765
    22182766!======================================================================
    22192767        SUBROUTINE interp_toga_time(day,day1,annee_ref
    … …  
    24863034        return
    24873035        end
     3036!======================================================================
     3037      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,
     3038     .       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
     3039!======================================================================
     3040      implicit none
     3041
     3042        integer nlev_max,kmax,kmax2
     3043        logical :: llesread = .true.
     3044
     3045        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
     3046     .  thlprof(nlev_max),
     3047     .  qprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
     3048     .  wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
     3049
     3050        integer, parameter :: ilesfile=1
     3051        integer :: ierr,irad,imax,jtot,k
     3052        logical :: lmoist,lcoriol,ltimedep
     3053        real :: xsize,ysize
     3054        real :: ustin,wsvsurf,timerad
     3055        character(80) :: chmess
     3056
     3057        if(.not.(llesread)) return
     3058
     3059       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
     3060        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
     3061        read (ilesfile,*) kmax
     3062        do k=1,kmax
     3063          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
     3064     .                      qprof (k),uprof(k),  vprof(k),  wprof(k),
     3065     .                      omega (k),o3mmr(k)
     3066        enddo
     3067        close(ilesfile)
     3068
     3069        return
     3070        end
     3071
     3072!======================================================================
     3073      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,
     3074     .    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
     3075!======================================================================
     3076      implicit none
     3077
     3078        integer nlev_max,kmax,kmax2
     3079        logical :: llesread = .true.
     3080
     3081        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
     3082     .  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),
     3083     .  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
     3084     .  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
     3085
     3086        integer, parameter :: ilesfile=1
     3087        integer :: ierr,irad,imax,jtot,k
     3088        logical :: lmoist,lcoriol,ltimedep
     3089        real :: xsize,ysize
     3090        real :: ustin,wsvsurf,timerad
     3091        character(80) :: chmess
     3092
     3093        if(.not.(llesread)) return
     3094
     3095       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
     3096        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
     3097        read (ilesfile,*) kmax
     3098        do k=1,kmax
     3099          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
     3100     .                qvprof (k),qlprof (k),qtprof (k),
     3101     .                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
     3102        enddo
     3103        close(ilesfile)
     3104
     3105        return
     3106        end
     3107
    24883108
    24893109
    … …  
    25873207      return
    25883208      end
    2589 
     3209!=====================================================================
     3210      subroutine read_amma(fich_amma,nlevel,ntime
     3211     :     ,zz,pp,temp,qv,u,v,dw
     3212     :     ,dt,dq,sens,flat)
     3213
     3214!program reading forcings of the AMMA case study
     3215
     3216
     3217      implicit none
     3218
     3219#include "netcdf.inc"
     3220
     3221      integer ntime,nlevel
     3222      integer l,k
     3223      character*80 :: fich_amma
     3224      real*8 time(ntime)
     3225      real*8 zz(nlevel)
     3226
     3227      real*8 temp(nlevel),pp(nlevel)
     3228      real*8 qv(nlevel),u(nlevel)
     3229      real*8 v(nlevel)
     3230      real*8 dw(nlevel,ntime)
     3231      real*8 dt(nlevel,ntime)
     3232      real*8 dq(nlevel,ntime)   
     3233      real*8 flat(ntime),sens(ntime)
     3234
     3235      integer nid, ierr
     3236      integer nbvar3d
     3237      parameter(nbvar3d=30)
     3238      integer var3didin(nbvar3d)
     3239
     3240      ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
     3241      if (ierr.NE.NF_NOERR) then
     3242         write(*,*) 'ERROR: Pb opening forcings nc file '
     3243         write(*,*) NF_STRERROR(ierr)
     3244         stop ""
     3245      endif
     3246
     3247
     3248       ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
     3249         if(ierr/=NF_NOERR) then
     3250           write(*,*) NF_STRERROR(ierr)
     3251           stop 'lev'
     3252         endif
     3253
     3254
     3255      ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
     3256         if(ierr/=NF_NOERR) then
     3257           write(*,*) NF_STRERROR(ierr)
     3258           stop 'temp'
     3259         endif
     3260
     3261      ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
     3262         if(ierr/=NF_NOERR) then
     3263           write(*,*) NF_STRERROR(ierr)
     3264           stop 'qv'
     3265         endif
     3266
     3267      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
     3268         if(ierr/=NF_NOERR) then
     3269           write(*,*) NF_STRERROR(ierr)
     3270           stop 'u'
     3271         endif
     3272
     3273      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
     3274         if(ierr/=NF_NOERR) then
     3275           write(*,*) NF_STRERROR(ierr)
     3276           stop 'v'
     3277         endif
     3278
     3279      ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
     3280         if(ierr/=NF_NOERR) then
     3281           write(*,*) NF_STRERROR(ierr)
     3282           stop 'dw'
     3283         endif
     3284
     3285      ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
     3286         if(ierr/=NF_NOERR) then
     3287           write(*,*) NF_STRERROR(ierr)
     3288           stop 'dt'
     3289         endif
     3290
     3291      ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
     3292         if(ierr/=NF_NOERR) then
     3293           write(*,*) NF_STRERROR(ierr)
     3294           stop 'dq'
     3295         endif
     3296     
     3297      ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
     3298         if(ierr/=NF_NOERR) then
     3299           write(*,*) NF_STRERROR(ierr)
     3300           stop 'sens'
     3301         endif
     3302
     3303      ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
     3304         if(ierr/=NF_NOERR) then
     3305           write(*,*) NF_STRERROR(ierr)
     3306           stop 'flat'
     3307         endif
     3308
     3309      ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
     3310         if(ierr/=NF_NOERR) then
     3311           write(*,*) NF_STRERROR(ierr)
     3312           stop 'pp'
     3313      endif
     3314
     3315!dimensions lecture
     3316!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
     3317 
     3318#ifdef NC_DOUBLE
     3319         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
     3320#else
     3321         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
     3322#endif
     3323         if(ierr/=NF_NOERR) then
     3324            write(*,*) NF_STRERROR(ierr)
     3325            stop "getvarup"
     3326         endif
     3327!          write(*,*)'lecture z ok',zz
     3328
     3329#ifdef NC_DOUBLE
     3330         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
     3331#else
     3332         ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
     3333#endif
     3334         if(ierr/=NF_NOERR) then
     3335            write(*,*) NF_STRERROR(ierr)
     3336            stop "getvarup"
     3337         endif
     3338!          write(*,*)'lecture th ok',temp
     3339
     3340#ifdef NC_DOUBLE
     3341         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
     3342#else
     3343         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
     3344#endif
     3345         if(ierr/=NF_NOERR) then
     3346            write(*,*) NF_STRERROR(ierr)
     3347            stop "getvarup"
     3348         endif
     3349!          write(*,*)'lecture qv ok',qv
     3350 
     3351#ifdef NC_DOUBLE
     3352         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
     3353#else
     3354         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
     3355#endif
     3356         if(ierr/=NF_NOERR) then
     3357            write(*,*) NF_STRERROR(ierr)
     3358            stop "getvarup"
     3359         endif
     3360!          write(*,*)'lecture u ok',u
     3361
     3362#ifdef NC_DOUBLE
     3363         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
     3364#else
     3365         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
     3366#endif
     3367         if(ierr/=NF_NOERR) then
     3368            write(*,*) NF_STRERROR(ierr)
     3369            stop "getvarup"
     3370         endif
     3371!          write(*,*)'lecture v ok',v
     3372
     3373#ifdef NC_DOUBLE
     3374         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
     3375#else
     3376         ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
     3377#endif
     3378         if(ierr/=NF_NOERR) then
     3379            write(*,*) NF_STRERROR(ierr)
     3380            stop "getvarup"
     3381         endif
     3382!          write(*,*)'lecture w ok',dw
     3383
     3384#ifdef NC_DOUBLE
     3385         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
     3386#else
     3387         ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
     3388#endif
     3389         if(ierr/=NF_NOERR) then
     3390            write(*,*) NF_STRERROR(ierr)
     3391            stop "getvarup"
     3392         endif
     3393!          write(*,*)'lecture dt ok',dt
     3394
     3395#ifdef NC_DOUBLE
     3396         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
     3397#else
     3398         ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
     3399#endif
     3400         if(ierr/=NF_NOERR) then
     3401            write(*,*) NF_STRERROR(ierr)
     3402            stop "getvarup"
     3403         endif
     3404!          write(*,*)'lecture dq ok',dq
     3405
     3406#ifdef NC_DOUBLE
     3407         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
     3408#else
     3409         ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
     3410#endif
     3411         if(ierr/=NF_NOERR) then
     3412            write(*,*) NF_STRERROR(ierr)
     3413            stop "getvarup"
     3414         endif
     3415!          write(*,*)'lecture sens ok',sens
     3416
     3417#ifdef NC_DOUBLE
     3418         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
     3419#else
     3420         ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
     3421#endif
     3422         if(ierr/=NF_NOERR) then
     3423            write(*,*) NF_STRERROR(ierr)
     3424            stop "getvarup"
     3425         endif
     3426!          write(*,*)'lecture flat ok',flat
     3427
     3428#ifdef NC_DOUBLE
     3429         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
     3430#else
     3431         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
     3432#endif
     3433         if(ierr/=NF_NOERR) then
     3434            write(*,*) NF_STRERROR(ierr)
     3435            stop "getvarup"
     3436         endif
     3437!          write(*,*)'lecture pp ok',pp
     3438
     3439         return
     3440         end subroutine read_amma
     3441!======================================================================
     3442        SUBROUTINE interp_amma_time(day,day1,annee_ref
     3443     i         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma
     3444     i         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
     3445     o         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
     3446        implicit none
     3447
     3448!---------------------------------------------------------------------------------------
     3449! Time interpolation of a 2D field to the timestep corresponding to day
     3450!
     3451! day: current julian day (e.g. 717538.2)
     3452! day1: first day of the simulation
     3453! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
     3454! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
     3455!---------------------------------------------------------------------------------------
     3456
     3457#include "compar1d.h"
     3458
     3459! inputs:
     3460        integer annee_ref
     3461        integer nt_amma,nlev_amma
     3462        integer year_ini_amma
     3463        real day, day1,day_ini_amma,dt_amma
     3464        real vitw_amma(nlev_amma,nt_amma)
     3465        real ht_amma(nlev_amma,nt_amma)
     3466        real hq_amma(nlev_amma,nt_amma)
     3467        real lat_amma(nt_amma)
     3468        real sens_amma(nt_amma)
     3469! outputs:
     3470        real vitw_prof(nlev_amma)
     3471        real ht_prof(nlev_amma)
     3472        real hq_prof(nlev_amma)
     3473        real lat_prof,sens_prof
     3474! local:
     3475        integer it_amma1, it_amma2,k
     3476        real timeit,time_amma1,time_amma2,frac
     3477
     3478
     3479        if (forcing_type.eq.6) then
     3480! Check that initial day of the simulation consistent with AMMA case:
     3481       if (annee_ref.ne.2006) then
     3482        print*,'Pour AMMA, annee_ref doit etre 2006'
     3483        print*,'Changer annee_ref dans run.def'
     3484        stop
     3485       endif
     3486       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
     3487        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
     3488        print*,'Changer dayref dans run.def'
     3489        stop
     3490       endif
     3491       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
     3492        print*,'AMMA a fini le 11 juillet'
     3493        print*,'Changer dayref ou nday dans run.def'
     3494        stop
     3495       endif
     3496       endif
     3497
     3498! Determine timestep relative to the 1st day of AMMA:
     3499!       timeit=(day-day1)*86400.
     3500!       if (annee_ref.eq.1992) then
     3501!        timeit=(day-day_ini_toga)*86400.
     3502!       else
     3503!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
     3504!       endif
     3505      timeit=(day-day_ini_amma)*86400
     3506
     3507! Determine the closest observation times:
     3508!       it_amma1=INT(timeit/dt_amma)+1
     3509!       it_amma2=it_amma1 + 1
     3510!       time_amma1=(it_amma1-1)*dt_amma
     3511!       time_amma2=(it_amma2-1)*dt_amma
     3512
     3513       it_amma1=INT(timeit/dt_amma)+1
     3514       IF (it_amma1 .EQ. nt_amma) THEN
     3515       it_amma2=it_amma1
     3516       ELSE
     3517       it_amma2=it_amma1 + 1
     3518       ENDIF
     3519       time_amma1=(it_amma1-1)*dt_amma
     3520       time_amma2=(it_amma2-1)*dt_amma
     3521
     3522       if (it_amma1 .gt. nt_amma) then
     3523        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '
     3524     :        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
     3525        stop
     3526       endif
     3527
     3528! time interpolation:
     3529       frac=(time_amma2-timeit)/(time_amma2-time_amma1)
     3530       frac=max(frac,0.0)
     3531
     3532       lat_prof = lat_amma(it_amma2)
     3533     :          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
     3534       sens_prof = sens_amma(it_amma2)
     3535     :          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
     3536
     3537       do k=1,nlev_amma
     3538        vitw_prof(k) = vitw_amma(k,it_amma2)
     3539     :          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
     3540        ht_prof(k) = ht_amma(k,it_amma2)
     3541     :          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
     3542        hq_prof(k) = hq_amma(k,it_amma2)
     3543     :          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
     3544        enddo
     3545
     3546        return
     3547        END
     3548
     3549!=====================================================================
     3550      subroutine read_fire(fich_fire,nlevel,ntime
     3551     :     ,zz,thl,qt,u,v,tke
     3552     :     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
     3553
     3554!program reading forcings of the FIRE case study
     3555
     3556
     3557      implicit none
     3558
     3559#include "netcdf.inc"
     3560
     3561      integer ntime,nlevel
     3562      integer l,k
     3563      character*80 :: fich_fire
     3564      real*8 time(ntime)
     3565      real*8 zz(nlevel)
     3566
     3567      real*8 thl(nlevel)
     3568      real*8 qt(nlevel),u(nlevel)
     3569      real*8 v(nlevel),tke(nlevel)
     3570      real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
     3571      real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
     3572      real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime) 
     3573
     3574      integer nid, ierr
     3575      integer nbvar3d
     3576      parameter(nbvar3d=30)
     3577      integer var3didin(nbvar3d)
     3578
     3579      ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)
     3580      if (ierr.NE.NF_NOERR) then
     3581         write(*,*) 'ERROR: Pb opening forcings nc file '
     3582         write(*,*) NF_STRERROR(ierr)
     3583         stop ""
     3584      endif
     3585
     3586
     3587       ierr=NF_INQ_VARID(nid,"zz",var3didin(1))
     3588         if(ierr/=NF_NOERR) then
     3589           write(*,*) NF_STRERROR(ierr)
     3590           stop 'lev'
     3591         endif
     3592
     3593
     3594      ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))
     3595         if(ierr/=NF_NOERR) then
     3596           write(*,*) NF_STRERROR(ierr)
     3597           stop 'temp'
     3598         endif
     3599
     3600      ierr=NF_INQ_VARID(nid,"qt",var3didin(3))
     3601         if(ierr/=NF_NOERR) then
     3602           write(*,*) NF_STRERROR(ierr)
     3603           stop 'qv'
     3604         endif
     3605
     3606      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
     3607         if(ierr/=NF_NOERR) then
     3608           write(*,*) NF_STRERROR(ierr)
     3609           stop 'u'
     3610         endif
     3611
     3612      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
     3613         if(ierr/=NF_NOERR) then
     3614           write(*,*) NF_STRERROR(ierr)
     3615           stop 'v'
     3616         endif
     3617
     3618      ierr=NF_INQ_VARID(nid,"tke",var3didin(6))
     3619         if(ierr/=NF_NOERR) then
     3620           write(*,*) NF_STRERROR(ierr)
     3621           stop 'tke'
     3622         endif
     3623
     3624      ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
     3625         if(ierr/=NF_NOERR) then
     3626           write(*,*) NF_STRERROR(ierr)
     3627           stop 'ug'
     3628         endif
     3629
     3630      ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
     3631         if(ierr/=NF_NOERR) then
     3632           write(*,*) NF_STRERROR(ierr)
     3633           stop 'vg'
     3634         endif
     3635     
     3636      ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
     3637         if(ierr/=NF_NOERR) then
     3638           write(*,*) NF_STRERROR(ierr)
     3639           stop 'wls'
     3640         endif
     3641
     3642      ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
     3643         if(ierr/=NF_NOERR) then
     3644           write(*,*) NF_STRERROR(ierr)
     3645           stop 'dqtdx'
     3646         endif
     3647
     3648      ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
     3649         if(ierr/=NF_NOERR) then
     3650           write(*,*) NF_STRERROR(ierr)
     3651           stop 'dqtdy'
     3652      endif
     3653
     3654      ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
     3655         if(ierr/=NF_NOERR) then
     3656           write(*,*) NF_STRERROR(ierr)
     3657           stop 'dqtdt'
     3658      endif
     3659
     3660      ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
     3661         if(ierr/=NF_NOERR) then
     3662           write(*,*) NF_STRERROR(ierr)
     3663           stop 'thl_rad'
     3664      endif
     3665!dimensions lecture
     3666!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
     3667 
     3668#ifdef NC_DOUBLE
     3669         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
     3670#else
     3671         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
     3672#endif
     3673         if(ierr/=NF_NOERR) then
     3674            write(*,*) NF_STRERROR(ierr)
     3675            stop "getvarup"
     3676         endif
     3677!          write(*,*)'lecture z ok',zz
     3678
     3679#ifdef NC_DOUBLE
     3680         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)
     3681#else
     3682         ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
     3683#endif
     3684         if(ierr/=NF_NOERR) then
     3685            write(*,*) NF_STRERROR(ierr)
     3686            stop "getvarup"
     3687         endif
     3688!          write(*,*)'lecture thl ok',thl
     3689
     3690#ifdef NC_DOUBLE
     3691         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
     3692#else
     3693         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
     3694#endif
     3695         if(ierr/=NF_NOERR) then
     3696            write(*,*) NF_STRERROR(ierr)
     3697            stop "getvarup"
     3698         endif
     3699!          write(*,*)'lecture qt ok',qt
     3700 
     3701#ifdef NC_DOUBLE
     3702         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
     3703#else
     3704         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
     3705#endif
     3706         if(ierr/=NF_NOERR) then
     3707            write(*,*) NF_STRERROR(ierr)
     3708            stop "getvarup"
     3709         endif
     3710!          write(*,*)'lecture u ok',u
     3711
     3712#ifdef NC_DOUBLE
     3713         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
     3714#else
     3715         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
     3716#endif
     3717         if(ierr/=NF_NOERR) then
     3718            write(*,*) NF_STRERROR(ierr)
     3719            stop "getvarup"
     3720         endif
     3721!          write(*,*)'lecture v ok',v
     3722
     3723#ifdef NC_DOUBLE
     3724         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)
     3725#else
     3726         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
     3727#endif
     3728         if(ierr/=NF_NOERR) then
     3729            write(*,*) NF_STRERROR(ierr)
     3730            stop "getvarup"
     3731         endif
     3732!          write(*,*)'lecture tke ok',tke
     3733
     3734#ifdef NC_DOUBLE
     3735         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
     3736#else
     3737         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)
     3738#endif
     3739         if(ierr/=NF_NOERR) then
     3740            write(*,*) NF_STRERROR(ierr)
     3741            stop "getvarup"
     3742         endif
     3743!          write(*,*)'lecture ug ok',ug
     3744
     3745#ifdef NC_DOUBLE
     3746         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)
     3747#else
     3748         ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)
     3749#endif
     3750         if(ierr/=NF_NOERR) then
     3751            write(*,*) NF_STRERROR(ierr)
     3752            stop "getvarup"
     3753         endif
     3754!          write(*,*)'lecture vg ok',vg
     3755
     3756#ifdef NC_DOUBLE
     3757         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)
     3758#else
     3759         ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
     3760#endif
     3761         if(ierr/=NF_NOERR) then
     3762            write(*,*) NF_STRERROR(ierr)
     3763            stop "getvarup"
     3764         endif
     3765!          write(*,*)'lecture wls ok',wls
     3766
     3767#ifdef NC_DOUBLE
     3768         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
     3769#else
     3770         ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
     3771#endif
     3772         if(ierr/=NF_NOERR) then
     3773            write(*,*) NF_STRERROR(ierr)
     3774            stop "getvarup"
     3775         endif
     3776!          write(*,*)'lecture dqtdx ok',dqtdx
     3777
     3778#ifdef NC_DOUBLE
     3779         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
     3780#else
     3781         ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
     3782#endif
     3783         if(ierr/=NF_NOERR) then
     3784            write(*,*) NF_STRERROR(ierr)
     3785            stop "getvarup"
     3786         endif
     3787!          write(*,*)'lecture dqtdy ok',dqtdy
     3788
     3789#ifdef NC_DOUBLE
     3790         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
     3791#else
     3792         ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
     3793#endif
     3794         if(ierr/=NF_NOERR) then
     3795            write(*,*) NF_STRERROR(ierr)
     3796            stop "getvarup"
     3797         endif
     3798!          write(*,*)'lecture dqtdt ok',dqtdt
     3799
     3800#ifdef NC_DOUBLE
     3801         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
     3802#else
     3803         ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
     3804#endif
     3805         if(ierr/=NF_NOERR) then
     3806            write(*,*) NF_STRERROR(ierr)
     3807            stop "getvarup"
     3808         endif
     3809!          write(*,*)'lecture thl_rad ok',thl_rad
     3810
     3811         return
     3812         end subroutine read_fire
Note: See TracChangeset for help on using the changeset viewer.