Changeset 1780 for LMDZ5/trunk/libf
- Timestamp:
- Jul 5, 2013, 10:38:13 AM (11 years ago)
- Location:
- LMDZ5/trunk/libf/phy1d
- Files:
-
- 2 added
- 5 edited
-
1DUTILS.h_no_writelim (modified) (11 diffs)
-
1D_decl_cases.h (modified) (2 diffs)
-
1D_interp_cases.h (modified) (2 diffs)
-
1D_read_forc_cases.h (modified) (8 diffs)
-
lmdz1d.F (modified) (17 diffs)
-
pbl_surface_mod.F90 (added)
-
tsoilnudge.h (added)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phy1d/1DUTILS.h_no_writelim
r1763 r1780 27 27 #include "compar1d.h" 28 28 #include "flux_arp.h" 29 #include "tsoilnudge.h" 29 30 #include "fcg_gcssold.h" 30 31 #include "fcg_racmo.h" … … 100 101 ! initial profiles from RICO idealized 101 102 ! LS convergence imposed from RICO (cst) 103 ! = 6 ==> forcing_amma = .true. 102 104 ! = 40 ==> forcing_GCSSold = .true. 103 105 ! initial profile from GCSS file 104 106 ! 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 105 121 ! 106 122 forcing_type = 0 … … 234 250 zpicinp = 300. 235 251 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 236 285 237 286 … … 257 306 write(lunout,*)' qsolinp = ', qsolinp 258 307 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 259 312 IF (forcing_type .eq.40) THEN 260 313 write(lunout,*) '--- Forcing type GCSS Old --- with:' … … 965 1018 SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 966 1019 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 ! 967 1024 ! Auteur : P. Le Van . 968 1025 ! … … 1409 1466 end 1410 1467 1468 c------------------------------------------------------------------------- 1469 SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) 1470 implicit none 1471 1472 c------------------------------------------------------------------------- 1473 c Read I.SANDU case forcing data 1474 c------------------------------------------------------------------------- 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 !===================================================================== 1508 c------------------------------------------------------------------------- 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 1513 c------------------------------------------------------------------------- 1514 c Read Astex case forcing data 1515 c------------------------------------------------------------------------- 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 1411 1552 !===================================================================== 1412 1553 subroutine read_twpice(fich_twpice,nlevel,ntime … … 1891 2032 return 1892 2033 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 2046 c------------------------------------------------------------------------- 2047 c Vertical interpolation of SANDUREF forcing data onto model levels 2048 c------------------------------------------------------------------------- 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 2128 cjyg 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 2163 c------------------------------------------------------------------------- 2164 c Vertical interpolation of Astex forcing data onto model levels 2165 c------------------------------------------------------------------------- 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 2249 cjyg 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 1893 2275 1894 2276 !====================================================================== … … 2055 2437 end 2056 2438 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 2057 2514 !===================================================================== 2058 2515 c------------------------------------------------------------------------- … … 2216 2673 end 2217 2674 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 2218 2766 !====================================================================== 2219 2767 SUBROUTINE interp_toga_time(day,day1,annee_ref … … 2486 3034 return 2487 3035 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 2488 3108 2489 3109 … … 2587 3207 return 2588 3208 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 -
LMDZ5/trunk/libf/phy1d/1D_decl_cases.h
r1607 r1780 81 81 real ht_proftwp(nlev_twpi),vt_proftwp(nlev_twpi) 82 82 real hq_proftwp(nlev_twpi),vq_proftwp(nlev_twpi) 83 84 85 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 86 !Declarations specifiques au cas AMMA 87 character*80 :: fich_amma 88 ! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp) 89 logical :: fixe_disvert=.true. 90 integer nlev_amma, nt_amma 91 ! parameter (nlev_amma=29, nt_amma=48) ! Fleur, juillet 2012 92 parameter (nlev_amma=36, nt_amma=48) ! Romain, octobre 2012 93 ! parameter (nlev_amma=26, nt_amma=48) ! Test MPL feverier 2013 94 integer year_ini_amma, day_ini_amma, mth_ini_amma 95 real heure_ini_amma 96 real day_ju_ini_amma ! Julian day of amma first day 97 parameter (year_ini_amma=2006) 98 parameter (mth_ini_amma=7) 99 parameter (day_ini_amma=10) ! 10 = 10Juil2006 100 parameter (heure_ini_amma=0.) !0h en secondes 101 real dt_amma 102 parameter (dt_amma=1800.) 103 104 !profils initiaux: 105 real plev_amma(nlev_amma) 106 real tv_amma(nlev_amma),rho_amma(nlev_amma) 107 real thv_amma(nlev_amma) 108 109 real z_amma(nlev_amma) 110 real th_amma(nlev_amma),q_amma(nlev_amma) 111 real u_amma(nlev_amma) 112 real v_amma(nlev_amma) 113 114 real thvsurf_amma,tvsurf_amma,rhosurf_amma,thsurf 115 116 real th_ammai(nlev_amma),q_ammai(nlev_amma) 117 real u_ammai(nlev_amma) 118 real v_ammai(nlev_amma) 119 real vitw_ammai(nlev_amma) 120 real ht_ammai(nlev_amma) 121 real hq_ammai(nlev_amma) 122 real vt_ammai(nlev_amma) 123 real vq_ammai(nlev_amma) 124 125 !forcings 126 real ht_amma(nlev_amma,nt_amma) 127 real hq_amma(nlev_amma,nt_amma) 128 real vitw_amma(nlev_amma,nt_amma) 129 real lat_amma(nt_amma),sens_amma(nt_amma) 130 131 !champs interpoles 132 real plev_profamma(nlev_amma),vitw_profamma(nlev_amma) 133 real ht_profamma(nlev_amma) 134 real hq_profamma(nlev_amma) 135 real lat_profamma,sens_profamma 136 real vt_profamma(nlev_amma) 137 real vq_profamma(nlev_amma) 138 real th_profamma(nlev_amma) 139 real q_profamma(nlev_amma) 140 real u_profamma(nlev_amma) 141 real v_profamma(nlev_amma) 142 143 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 144 !Declarations specifiques au cas FIRE 145 character*80 :: fich_fire 146 integer nlev_fire, nt_fire 147 parameter (nlev_fire=120, nt_fire=1) 148 integer year_ini_fire, day_ini_fire, mth_ini_fire 149 real heure_ini_fire 150 real day_ju_ini_fire ! Julian day of fire first day 151 parameter (year_ini_fire=1987) 152 parameter (mth_ini_fire=7) 153 parameter (day_ini_fire=14) ! 14 = 14Juil1987 154 parameter (heure_ini_fire=0.) !0h en secondes 155 156 !profils initiaux: 157 real z_fire(nlev_fire) 158 real thl_fire(nlev_fire),qt_fire(nlev_fire) 159 real u_fire(nlev_fire), v_fire(nlev_fire) 160 real tke_fire(nlev_fire) 161 162 !forcings 163 real ugeo_fire(nlev_fire),vgeo_fire(nlev_fire) 164 real wls_fire(nlev_fire),dqtdx_fire(nlev_fire) 165 real dqtdy_fire(nlev_fire) 166 real dqtdt_fire(nlev_fire),thl_rad_fire(nlev_fire) 167 168 real ugeo_mod(llm),vgeo_mod(llm),wls_mod(llm) 169 real dqtdx_mod(llm),dqtdy_mod(llm),dqtdt_mod(llm) 170 real thl_rad_mod(llm) 83 171 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 84 172 ! Declarations specifiques au cas GCSSold … … 127 215 real sens_prof,flat_prof,fact 128 216 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 129 217 ! declarations specifiques au cas Sandu 218 character*80 :: fich_sandu 219 ! integer nlev_prof 220 ! parameter (nlev_prof = 41) 221 integer nlev_sandu, nt_sandu 222 parameter (nlev_sandu=87, nt_sandu=13) 223 integer year_ini_sandu, day_ini_sandu, mth_ini_sandu 224 real day_ju_ini_sandu ! Julian day of sandu case first day 225 parameter (year_ini_sandu=2006) 226 parameter (mth_ini_sandu=7) 227 parameter (day_ini_sandu=15) ! 196 = 15 juillet 2006 228 real dt_sandu, tau_sandu 229 logical :: trouve_700=.true. 230 parameter (dt_sandu=6.*3600.) ! forcages donnes ttes les 6 heures par ifa_sandu.txt 231 ! parameter (tau_sandu=3600.) ! temps de relaxation u,v,thetal,qt vers profil init et au dessus 700hPa 232 !! 233 integer it_sandu1, it_sandu2 234 real time_sandu1,time_sandu2 235 236 real ts_sandu(nt_sandu) 237 ! profs comme "profil sandu" 238 real plev_profs(nlev_sandu) 239 real t_profs(nlev_sandu),thl_profs(nlev_sandu) 240 real q_profs(nlev_sandu) 241 real u_profs(nlev_sandu),v_profs(nlev_sandu),w_profs(nlev_sandu) 242 real omega_profs(nlev_sandu),o3mmr_profs(nlev_sandu) 243 244 real thl_mod(llm),omega_mod(llm),o3mmr_mod(llm),tke_mod(llm) 245 ! pour relaxer u,v,thl et qt vers les profils initiaux au dessus de 700hPa 246 real relax_u(llm),relax_v(llm),relax_thl(llm),relax_q(llm,2) 247 !vertical advection computation 248 real d_t_z(llm), d_q_z(llm) 249 real d_t_dyn_z(llm), d_q_dyn_z(llm) 250 real zz(llm) 251 real zfact 252 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 253 ! Declarations specifiques au cas Astex 254 character*80 :: fich_astex 255 integer nlev_astex, nt_astex 256 parameter (nlev_astex=34, nt_astex=49) 257 integer year_ini_astex, day_ini_astex, mth_ini_astex 258 real day_ju_ini_astex ! Julian day of astex case first day 259 parameter (year_ini_astex=1992) 260 parameter (mth_ini_astex=6) 261 parameter (day_ini_astex=13) ! 165 = 13 juin 1992 262 real dt_astex, tau_astex 263 parameter (dt_astex=3600.) ! forcages donnes ttes les heures par ifa_astex.txt 264 integer it_astex1, it_astex2 265 real time_astex1,time_astex2 266 real ts_astex(nt_astex),div_astex(nt_astex),ug_astex(nt_astex) 267 real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex) 268 real div_prof,ug_prof,vg_prof,ufa_prof,vfa_prof 269 ! profa comme "profil astex" 270 real plev_profa(nlev_astex) 271 real t_profa(nlev_astex),thl_profa(nlev_astex) 272 real qv_profa(nlev_astex),ql_profa(nlev_astex) 273 real qt_profa(nlev_astex),o3mmr_profa(nlev_astex) 274 real u_profa(nlev_astex),v_profa(nlev_astex),w_profa(nlev_astex) 275 real tke_profa(nlev_astex) 276 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 277 -
LMDZ5/trunk/libf/phy1d/1D_interp_cases.h
r1607 r1780 181 181 182 182 endif ! forcing_twpice 183 184 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 185 !--------------------------------------------------------------------- 186 ! Interpolation forcing AMMA 187 !--------------------------------------------------------------------- 188 189 if (forcing_amma) then 190 191 print*, 192 : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=', 193 : daytime,day1,(daytime-day1)*86400., 194 : (daytime-day1)*86400/dt_amma 195 196 ! time interpolation using TOGA interpolation routine 197 CALL interp_amma_time(daytime,day1,annee_ref 198 i ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma 199 i ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma 200 o ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma 201 : ,sens_profamma) 202 203 print*,'apres interpolation temporelle AMMA' 204 205 do k=1,nlev_amma 206 th_profamma(k)=0. 207 q_profamma(k)=0. 208 u_profamma(k)=0. 209 v_profamma(k)=0. 210 vt_profamma(k)=0. 211 vq_profamma(k)=0. 212 enddo 213 ! vertical interpolation using TOGA interpolation routine: 214 ! write(*,*)'avant interp vert', t_proftwp 215 CALL interp_toga_vertical(play,nlev_amma,plev_amma 216 : ,th_profamma,q_profamma,u_profamma,v_profamma 217 : ,vitw_profamma 218 : ,ht_profamma,vt_profamma,hq_profamma,vq_profamma 219 : ,t_mod,q_mod,u_mod,v_mod,w_mod 220 : ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) 221 write(*,*) 'Profil initial forcing AMMA interpole' 222 223 224 !calcul de l'advection verticale a partir du omega 225 cCalcul des gradients verticaux 226 cinitialisation 227 do l=1,llm 228 d_t_z(l)=0. 229 d_q_z(l)=0. 230 enddo 231 232 DO l=2,llm-1 233 d_t_z(l)=(temp(l+1)-temp(l-1)) 234 & /(play(l+1)-play(l-1)) 235 d_q_z(l)=(q(l+1,1)-q(l-1,1)) 236 & /(play(l+1)-play(l-1)) 237 ENDDO 238 d_t_z(1)=d_t_z(2) 239 d_q_z(1)=d_q_z(2) 240 d_t_z(llm)=d_t_z(llm-1) 241 d_q_z(llm)=d_q_z(llm-1) 242 243 244 do l = 1, llm 245 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 246 omega(l) = w_mod(l)*(-rg*rho(l)) 247 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 248 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 249 !calcul de l'advection totale 250 ! d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l) 251 !attention: on impose dth 252 d_th_adv(l) = alpha*omega(l)/rcpd+ 253 & ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l) 254 ! d_th_adv(l) = 0. 255 ! print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l) 256 d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l) 257 ! d_q_adv(l,1) = 0. 258 ! print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l) 259 260 dt_cooling(l) = 0.0 261 enddo 262 263 264 ! ok_flux_surf=.false. 265 fsens=-1.*sens_profamma 266 flat=-1.*lat_profamma 267 268 endif ! forcing_amma 269 183 270 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 184 271 !--------------------------------------------------------------------- … … 254 341 endif ! forcing_armcu 255 342 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 256 343 !--------------------------------------------------------------------- 344 ! Interpolation forcing in time and onto model levels 345 !--------------------------------------------------------------------- 346 if (forcing_sandu) then 347 348 print*, 349 : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=', 350 : day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu 351 352 ! time interpolation: 353 ! ATTENTION, cet appel ne convient pas pour TOGA !! 354 ! revoir 1DUTILS.h et les arguments 355 CALL interp_sandu_time(daytime,day1,annee_ref 356 i ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu 357 i ,nlev_sandu 358 i ,ts_sandu,ts_prof) 359 360 if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d 361 362 ! vertical interpolation: 363 CALL interp_sandu_vertical(play,nlev_sandu,plev_profs 364 : ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs 365 : ,omega_profs,o3mmr_profs 366 : ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod 367 : ,omega_mod,o3mmr_mod,mxcalc) 368 !calcul de l'advection verticale 369 cCalcul des gradients verticaux 370 cinitialisation 371 d_t_z(:)=0. 372 d_q_z(:)=0. 373 d_t_dyn_z(:)=0. 374 d_q_dyn_z(:)=0. 375 ! schema centre 376 ! DO l=2,llm-1 377 ! d_t_z(l)=(temp(l+1)-temp(l-1)) 378 ! & /(play(l+1)-play(l-1)) 379 ! d_q_z(l)=(q(l+1,1)-q(l-1,1)) 380 ! & /(play(l+1)-play(l-1)) 381 ! schema amont 382 DO l=2,llm-1 383 d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l)) 384 d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l)) 385 ! print *,'l temp2 temp0 play2 play0 omega_mod', 386 ! & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l) 387 ENDDO 388 d_t_z(1)=d_t_z(2) 389 d_q_z(1)=d_q_z(2) 390 d_t_z(llm)=d_t_z(llm-1) 391 d_q_z(llm)=d_q_z(llm-1) 392 393 ! calcul de l advection verticale 394 ! Confusion w (m/s) et omega (Pa/s) !! 395 d_t_dyn_z(:)=omega_mod(:)*d_t_z(:) 396 d_q_dyn_z(:)=omega_mod(:)*d_q_z(:) 397 ! do l=1,llm 398 ! print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z', 399 ! :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l) 400 ! enddo 401 402 403 ! large-scale forcing : pour le cas Sandu ces forcages sont la SST 404 ! et une divergence constante -> profil de omega 405 tsurf = ts_prof 406 write(*,*) 'SST suivante: ',tsurf 407 do l = 1, llm 408 omega(l) = omega_mod(l) 409 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 410 411 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 412 ! 413 ! d_th_adv(l) = 0.0 414 ! d_q_adv(l,1) = 0.0 415 !CR:test advection=0 416 !calcul de l'advection verticale 417 d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l) 418 ! print*,'temp adv',l,-d_t_dyn_z(l) 419 d_q_adv(l,1) = -d_q_dyn_z(l) 420 ! print*,'q adv',l,-d_q_dyn_z(l) 421 dt_cooling(l) = 0.0 422 enddo 423 endif ! forcing_sandu 424 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 425 !--------------------------------------------------------------------- 426 ! Interpolation forcing in time and onto model levels 427 !--------------------------------------------------------------------- 428 if (forcing_astex) then 429 430 print*, 431 : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=', 432 : day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex 433 434 ! time interpolation: 435 ! ATTENTION, cet appel ne convient pas pour TOGA !! 436 ! revoir 1DUTILS.h et les arguments 437 CALL interp_astex_time(daytime,day1,annee_ref 438 i ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex 439 i ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex 440 i ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof 441 i ,ufa_prof,vfa_prof) 442 443 if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d 444 445 ! vertical interpolation: 446 CALL interp_astex_vertical(play,nlev_astex,plev_profa 447 : ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa 448 : ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa 449 : ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod 450 : ,tke_mod,o3mmr_mod,mxcalc) 451 !calcul de l'advection verticale 452 !Calcul des gradients verticaux 453 !initialisation 454 d_t_z(:)=0. 455 d_q_z(:)=0. 456 d_t_dyn_z(:)=0. 457 d_q_dyn_z(:)=0. 458 ! schema centre 459 ! DO l=2,llm-1 460 ! d_t_z(l)=(temp(l+1)-temp(l-1)) 461 ! & /(play(l+1)-play(l-1)) 462 ! d_q_z(l)=(q(l+1,1)-q(l-1,1)) 463 ! & /(play(l+1)-play(l-1)) 464 ! schema amont 465 DO l=2,llm-1 466 d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l)) 467 d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l)) 468 ! print *,'l temp2 temp0 play2 play0 omega_mod', 469 ! & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l) 470 ENDDO 471 d_t_z(1)=d_t_z(2) 472 d_q_z(1)=d_q_z(2) 473 d_t_z(llm)=d_t_z(llm-1) 474 d_q_z(llm)=d_q_z(llm-1) 475 476 ! calcul de l advection verticale 477 ! Confusion w (m/s) et omega (Pa/s) !! 478 d_t_dyn_z(:)=w_mod(:)*d_t_z(:) 479 d_q_dyn_z(:)=w_mod(:)*d_q_z(:) 480 ! do l=1,llm 481 ! print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z', 482 ! :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l) 483 ! enddo 484 485 486 ! large-scale forcing : pour le cas Astex ces forcages sont la SST 487 ! la divergence,ug,vg,ufa,vfa 488 tsurf = ts_prof 489 write(*,*) 'SST suivante: ',tsurf 490 do l = 1, llm 491 omega(l) = w_mod(l) 492 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 493 494 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 495 ! 496 ! d_th_adv(l) = 0.0 497 ! d_q_adv(l,1) = 0.0 498 !CR:test advection=0 499 !calcul de l'advection verticale 500 d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l) 501 ! print*,'temp adv',l,-d_t_dyn_z(l) 502 d_q_adv(l,1) = -d_q_dyn_z(l) 503 ! print*,'q adv',l,-d_q_dyn_z(l) 504 dt_cooling(l) = 0.0 505 enddo 506 endif ! forcing_astex 507 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 508 -
LMDZ5/trunk/libf/phy1d/1D_read_forc_cases.h
r1607 r1780 4 4 !---------------------------------------------------------------------- 5 5 6 if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then 7 6 if (forcing_les .or. forcing_radconv 7 : .or. forcing_GCSSold .or. forcing_fire) then 8 9 if (forcing_fire) then 10 !---------------------------------------------------------------------- 11 !read fire forcings from fire.nc 12 !---------------------------------------------------------------------- 13 fich_fire='fire.nc' 14 call read_fire(fich_fire,nlev_fire,nt_fire 15 : ,height,tttprof,qtprof,uprof,vprof,e12prof 16 : ,ugprof,vgprof,wfls,dqtdxls 17 : ,dqtdyls,dqtdtls,thlpcar) 18 write(*,*) 'Forcing FIRE lu' 19 kmax=120 ! nombre de niveaux dans les profils et forcages 20 else 8 21 !---------------------------------------------------------------------- 9 22 ! Read profiles from files: prof.inp.001 and lscale.inp.001 … … 16 29 . wfls,dqtdxls,dqtdyls,dqtdtls, 17 30 . thlpcar) 31 endif 18 32 19 33 ! compute altitudes of play levels. … … 35 49 frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1)) 36 50 ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1)) 37 if (forcing_GCSSold .AND. tp_ini_GCSSold)then ! pot. temp. in initial profile51 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 38 52 temp(l) = ttt*(play(l)/pzero)**rkappa 39 53 teta(l) = ttt 40 54 else 41 55 temp(l) = ttt 42 56 teta(l) = ttt*(pzero/play(l))**rkappa 43 57 endif 44 58 print *,' temp,teta ',l,temp(l),teta(l) 45 59 q(l,1) = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1)) … … 58 72 if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then 59 73 ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1)) 60 if (forcing_GCSSold .AND. tp_ini_GCSSold)then ! pot. temp. in initial profile74 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 61 75 temp(l) = ttt*(play(l)/pzero)**rkappa 62 76 teta(l) = ttt 63 77 else 64 78 temp(l) = ttt 65 79 teta(l) = ttt*(pzero/play(l))**rkappa 66 80 endif 67 81 print *,' temp,teta ',l,temp(l),teta(l) 68 82 q(l,1) = qtprof(k)-frac*( qtprof(k)- qtprof(k-1)) … … 77 91 elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1) 78 92 ttt =tttprof(1) 79 if (forcing_GCSSold .AND. tp_ini_GCSSold)then ! pot. temp. in initial profile93 if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile 80 94 temp(l) = ttt*(play(l)/pzero)**rkappa 81 95 teta(l) = ttt 82 96 else 83 97 temp(l) = ttt 84 98 teta(l) = ttt*(pzero/play(l))**rkappa 85 99 endif 86 100 q(l,1) = qtprof(1) 87 101 u(l) = uprof(1) … … 112 126 enddo 113 127 114 endif ! forcing_les .or. forcing_GCSSold 128 endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire 115 129 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 116 130 !--------------------------------------------------------------------- … … 263 277 264 278 endif !forcing_twpice 279 280 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 281 !--------------------------------------------------------------------- 282 ! Forcing from AMMA experiment (Couvreux et al. 2010) : 283 !--------------------------------------------------------------------- 284 285 if (forcing_amma) then 286 !read AMMA forcings 287 fich_amma='amma.nc' 288 call read_amma(fich_amma,nlev_amma,nt_amma 289 : ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma 290 : ,ht_amma,hq_amma,sens_amma,lat_amma) 291 292 write(*,*) 'Forcing AMMA lu' 293 294 !champs initiaux: 295 do k=1,nlev_amma 296 th_ammai(k)=th_amma(k) 297 q_ammai(k)=q_amma(k) 298 u_ammai(k)=u_amma(k) 299 v_ammai(k)=v_amma(k) 300 vitw_ammai(k)=vitw_amma(k,12) 301 ht_ammai(k)=ht_amma(k,12) 302 hq_ammai(k)=hq_amma(k,12) 303 vt_ammai(k)=0. 304 vq_ammai(k)=0. 305 enddo 306 omega(:)=0. 307 omega2(:)=0. 308 rho(:)=0. 309 ! vertical interpolation using TOGA interpolation routine: 310 ! write(*,*)'avant interp vert', t_proftwp 311 CALL interp_toga_vertical(play,nlev_amma,plev_amma 312 : ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai 313 : ,ht_ammai,vt_ammai,hq_ammai,vq_ammai 314 : ,t_mod,q_mod,u_mod,v_mod,w_mod 315 : ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) 316 ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod 317 318 ! initial and boundary conditions : 319 ! tsurf = ts_proftwp 320 write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc 321 do l = 1, llm 322 ! Ligne du dessous à decommenter si on lit theta au lieu de temp 323 ! temp(l) = t_mod(l)*(play(l)/pzero)**rkappa 324 temp(l) = t_mod(l) 325 q(l,1) = q_mod(l) 326 q(l,2) = 0.0 327 ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) 328 u(l) = u_mod(l) 329 v(l) = v_mod(l) 330 rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 331 omega(l) = w_mod(l)*(-rg*rho(l)) 332 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 333 334 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 335 !on applique le forcage total au premier pas de temps 336 !attention: signe different de toga 337 d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 338 !forcage en th 339 ! d_th_adv(l) = ht_mod(l) 340 d_q_adv(l,1) = hq_mod(l) 341 d_q_adv(l,2) = 0.0 342 dt_cooling(l)=0. 343 enddo 344 write(*,*) 'Profil initial forcing AMMA interpole temp39', 345 & temp(39) 346 347 348 ! ok_flux_surf=.false. 349 fsens=-1.*sens_amma(12) 350 flat=-1.*lat_amma(12) 351 352 endif !forcing_amma 353 354 265 355 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 266 356 !--------------------------------------------------------------------- … … 366 456 endif ! forcing_armcu 367 457 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 368 458 !--------------------------------------------------------------------- 459 ! Forcing from transition case of Irina Sandu 460 !--------------------------------------------------------------------- 461 462 if (forcing_sandu) then 463 write(*,*) 'Avant lecture Forcing SANDU' 464 465 ! read sanduref forcing : 466 fich_sandu = './ifa_sanduref.txt' 467 CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) 468 469 write(*,*) 'Forcing SANDU lu' 470 471 !---------------------------------------------------------------------- 472 ! Read profiles from file: prof.inp.001 473 !---------------------------------------------------------------------- 474 475 call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs, 476 . thl_profs,q_profs,u_profs,v_profs, 477 . w_profs,omega_profs,o3mmr_profs) 478 479 ! time interpolation for initial conditions: 480 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 481 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! 482 ! revoir 1DUTILS.h et les arguments 483 484 print *,'Avant interp_sandu_time' 485 print *,'daytime=',daytime 486 print *,'day1=',day1 487 print *,'annee_ref=',annee_ref 488 print *,'year_ini_sandu=',year_ini_sandu 489 print *,'day_ju_ini_sandu=',day_ju_ini_sandu 490 print *,'nt_sandu=',nt_sandu 491 print *,'dt_sandu=',dt_sandu 492 print *,'nlev_sandu=',nlev_sandu 493 CALL interp_sandu_time(daytime,day1,annee_ref 494 i ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu 495 i ,nlev_sandu 496 i ,ts_sandu,ts_prof) 497 498 ! vertical interpolation: 499 print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu 500 CALL interp_sandu_vertical(play,nlev_sandu,plev_profs 501 : ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs 502 : ,omega_profs,o3mmr_profs 503 : ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod 504 : ,omega_mod,o3mmr_mod,mxcalc) 505 write(*,*) 'Profil initial forcing SANDU interpole' 506 507 ! initial and boundary conditions : 508 tsurf = ts_prof 509 write(*,*) 'SST initiale: ',tsurf 510 do l = 1, llm 511 temp(l) = t_mod(l) 512 tetal(l)=thl_mod(l) 513 q(l,1) = q_mod(l) 514 q(l,2) = 0.0 515 u(l) = u_mod(l) 516 v(l) = v_mod(l) 517 w(l) = w_mod(l) 518 omega(l) = omega_mod(l) 519 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 520 !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 521 !? omega2(l)=-rho(l)*omega(l) 522 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 523 ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 524 ! d_q_adv(l,1) = vq_mod(l) 525 d_th_adv(l) = alpha*omega(l)/rcpd 526 d_q_adv(l,1) = 0.0 527 d_q_adv(l,2) = 0.0 528 enddo 529 530 endif ! forcing_sandu 531 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 532 !--------------------------------------------------------------------- 533 ! Forcing from Astex case 534 !--------------------------------------------------------------------- 535 536 if (forcing_astex) then 537 write(*,*) 'Avant lecture Forcing Astex' 538 539 ! read astex forcing : 540 fich_astex = './ifa_astex.txt' 541 CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex, 542 : ug_astex,vg_astex,ufa_astex,vfa_astex) 543 544 write(*,*) 'Forcing Astex lu' 545 546 !---------------------------------------------------------------------- 547 ! Read profiles from file: prof.inp.001 548 !---------------------------------------------------------------------- 549 550 call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa, 551 . thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa, 552 . w_profa,tke_profa,o3mmr_profa) 553 554 ! time interpolation for initial conditions: 555 write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 556 ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! 557 ! revoir 1DUTILS.h et les arguments 558 559 print *,'Avant interp_astex_time' 560 print *,'daytime=',daytime 561 print *,'day1=',day1 562 print *,'annee_ref=',annee_ref 563 print *,'year_ini_astex=',year_ini_astex 564 print *,'day_ju_ini_astex=',day_ju_ini_astex 565 print *,'nt_astex=',nt_astex 566 print *,'dt_astex=',dt_astex 567 print *,'nlev_astex=',nlev_astex 568 CALL interp_astex_time(daytime,day1,annee_ref 569 i ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex 570 i ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex 571 i ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof 572 i ,ufa_prof,vfa_prof) 573 574 ! vertical interpolation: 575 print *,'Avant interp_vertical: nlev_astex=',nlev_astex 576 CALL interp_astex_vertical(play,nlev_astex,plev_profa 577 : ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa 578 : ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa 579 : ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod 580 : ,tke_mod,o3mmr_mod,mxcalc) 581 write(*,*) 'Profil initial forcing Astex interpole' 582 583 ! initial and boundary conditions : 584 tsurf = ts_prof 585 write(*,*) 'SST initiale: ',tsurf 586 do l = 1, llm 587 temp(l) = t_mod(l) 588 tetal(l)=thl_mod(l) 589 q(l,1) = qv_mod(l) 590 q(l,2) = ql_mod(l) 591 u(l) = u_mod(l) 592 v(l) = v_mod(l) 593 w(l) = w_mod(l) 594 omega(l) = w_mod(l) 595 ! omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 596 ! rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) 597 ! omega2(l)=-rho(l)*omega(l) 598 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 599 ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 600 ! d_q_adv(l,1) = vq_mod(l) 601 d_th_adv(l) = alpha*omega(l)/rcpd 602 d_q_adv(l,1) = 0.0 603 d_q_adv(l,2) = 0.0 604 enddo 605 606 endif ! forcing_astex 607 -
LMDZ5/trunk/libf/phy1d/lmdz1d.F
r1763 r1780 6 6 use dimphy 7 7 use surface_data, only : type_ocean,ok_veget 8 use pbl_surface_mod, only : pbl_surface_init, pbl_surface_final 8 use pbl_surface_mod, only : ftsoil, pbl_surface_init, 9 $ pbl_surface_final 9 10 use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 10 11 … … 25 26 #include "compar1d.h" 26 27 #include "flux_arp.h" 28 #include "tsoilnudge.h" 27 29 #include "fcg_gcssold.h" 28 30 !!!#include "fbforcing.h" … … 86 88 87 89 integer :: kmax = llm 88 integer nlev_max 89 parameter (nlev_max = 100 )90 integer nlev_max,llm700 91 parameter (nlev_max = 1000) 90 92 real timestep, frac, timeit 91 93 real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max), … … 98 100 c integer :: forcing_type 99 101 logical :: forcing_les = .false. 100 logical :: forcing_armcu = .false.102 logical :: forcing_armcu = .false. 101 103 logical :: forcing_rico = .false. 102 104 logical :: forcing_radconv = .false. 103 105 logical :: forcing_toga = .false. 104 106 logical :: forcing_twpice = .false. 107 logical :: forcing_amma = .false. 105 108 logical :: forcing_GCM2SCM = .false. 106 109 logical :: forcing_GCSSold = .false. 110 logical :: forcing_sandu = .false. 111 logical :: forcing_astex = .false. 112 logical :: forcing_fire = .false. 107 113 integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file 108 114 ! (cf read_tsurf1d.F) 109 115 110 116 !vertical advection computation 111 112 113 114 117 ! real d_t_z(llm), d_q_z(llm) 118 ! real d_t_dyn_z(llm), d_q_dyn_z(llm) 119 ! real zz(llm) 120 ! real zfact 115 121 116 122 !flag forcings … … 129 135 real :: pzero=1.e5 130 136 real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1) 131 real :: playd(llm),zlayd(llm) 137 real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub 132 138 133 139 !--------------------------------------------------------------------- … … 137 143 integer :: iq 138 144 real :: phi(llm) 145 real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm) 139 146 real :: rlat_rad(1),rlon_rad(1) 140 real :: teta(llm),temp(llm),u(llm),v(llm)141 147 real :: omega(llm+1),omega2(llm),rho(llm+1) 142 148 real :: ug(llm),vg(llm),fcoriolis … … 196 202 ! Fichiers et d'autres variables 197 203 !--------------------------------------------------------------------- 198 real ttt 204 real ttt,bow,q1 199 205 integer :: ierr,k,l,i,it=1,mxcalc 200 206 integer jjmp1 … … 252 258 ! initial profiles from RICO files 253 259 ! LS convergence imposed from RICO files 260 !forcing_type = 6 ==> forcing_amma = .true. 261 ! initial profiles from AMMA nc file 262 ! LS convergence, omega and surface fluxes imposed from AMMA file 254 263 !forcing_type = 40 ==> forcing_GCSSold = .true. 255 264 ! initial profile from GCSS file 256 265 ! LS convergence imposed from GCSS file 266 !forcing_type = 59 ==> forcing_sandu = .true. 267 ! initial profiles from sanduref file: see prof.inp.001 268 ! SST varying with time and divergence constante: see ifa_sanduref.txt file 269 ! Radiation has to be computed interactively 270 !forcing_type = 60 ==> forcing_astex = .true. 271 ! initial profiles from file: see prof.inp.001 272 ! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file 273 ! Radiation has to be computed interactively 257 274 !forcing_type = 61 ==> forcing_armcu = .true. 258 ! initial profile from arm_cu file 259 ! LS convergence imposed from arm_cu file 275 ! initial profiles from file: see prof.inp.001 276 ! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt 277 ! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt 278 ! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 279 ! Radiation to be switched off 260 280 ! 261 281 if (forcing_type .eq.0) THEN … … 271 291 elseif (forcing_type .eq.5) THEN 272 292 forcing_rico = .true. 293 elseif (forcing_type .eq.6) THEN 294 forcing_amma = .true. 273 295 elseif (forcing_type .eq.40) THEN 274 296 forcing_GCSSold = .true. 297 elseif (forcing_type .eq.59) THEN 298 forcing_sandu = .true. 299 elseif (forcing_type .eq.60) THEN 300 forcing_astex = .true. 275 301 elseif (forcing_type .eq.61) THEN 276 302 forcing_armcu = .true. … … 278 304 else 279 305 write (*,*) 'ERROR : unknown forcing_type ', forcing_type 280 stop 'Forcing_type should be 0,1,2,3 or 40'306 stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61' 281 307 ENDIF 282 308 print*,"forcing type=",forcing_type … … 288 314 289 315 type_ts_forcing = 0 290 if (forcing_toga) type_ts_forcing = 1 316 if (forcing_toga .or. forcing_sandu .or. forcing_astex) 317 : type_ts_forcing = 1 291 318 292 319 !--------------------------------------------------------------------- … … 327 354 c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026) 328 355 IF(forcing_type .EQ. 61) fnday=53100./86400. 356 c Special case for amma which lasts less than one day : 64800s !! (MPL 20120216) 357 IF(forcing_type .EQ. 6) fnday=64800./86400. 329 358 annee_ref = anneeref 330 359 mois = 1 … … 336 365 day_ini = day 337 366 day_end = day_ini + nday 367 368 IF (forcing_type .eq.2) THEN 338 369 ! Convert the initial date of Toga-Coare to Julian day 339 370 call ymds2ju 340 371 $ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga) 341 372 373 ELSEIF (forcing_type .eq.4) THEN 342 374 ! Convert the initial date of TWPICE to Julian day 343 375 call ymds2ju 344 376 $ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi 345 377 $ ,day_ju_ini_twpi) 346 347 ! Convert the initial date of Arm_cu to Julian day 378 ELSEIF (forcing_type .eq.6) THEN 379 ! Convert the initial date of AMMA to Julian day 380 call ymds2ju 381 $ (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma 382 $ ,day_ju_ini_amma) 383 384 ELSEIF (forcing_type .eq.59) THEN 385 ! Convert the initial date of Sandu case to Julian day 386 call ymds2ju 387 $ (year_ini_sandu,mth_ini_sandu,day_ini_sandu, 388 $ time_ini*3600.,day_ju_ini_sandu) 389 390 ELSEIF (forcing_type .eq.60) THEN 391 ! Convert the initial date of Astex case to Julian day 392 call ymds2ju 393 $ (year_ini_astex,mth_ini_astex,day_ini_astex, 394 $ time_ini*3600.,day_ju_ini_astex) 395 396 ELSEIF (forcing_type .eq.61) THEN 397 398 ! Convert the initial date of Arm_cu case to Julian day 348 399 call ymds2ju 349 400 $ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu 350 401 $ ,day_ju_ini_armcu) 402 ENDIF 351 403 352 404 daytime = day + time_ini/24. ! 1st day and initial time of the simulation … … 435 487 ccc zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles 436 488 489 IF (forcing_type .eq. 59) THEN 490 ! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m 437 491 write(*,*) '***********************' 438 492 do l = 1, llm 439 493 write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l) 494 if (trouve_700 .and. play(l).le.70000) then 495 llm700=l 496 print *,'llm700,play=',llm700,play(l)/100. 497 trouve_700= .false. 498 endif 440 499 enddo 441 500 write(*,*) '***********************' 501 ENDIF 442 502 443 503 c … … 515 575 agesno = xagesno 516 576 tsoil(:,:,:)=tsurf 577 !------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012) 578 ! tsoil(1,1,1)=299.18 579 ! tsoil(1,2,1)=300.08 580 ! tsoil(1,3,1)=301.88 581 ! tsoil(1,4,1)=305.48 582 ! tsoil(1,5,1)=308.00 583 ! tsoil(1,6,1)=308.00 584 ! tsoil(1,7,1)=308.00 585 ! tsoil(1,8,1)=308.00 586 ! tsoil(1,9,1)=308.00 587 ! tsoil(1,10,1)=308.00 588 ! tsoil(1,11,1)=308.00 589 !----------------------------------------------------------------------- 517 590 call pbl_surface_init(qsol, fder, snsrf, qsurfsrf, 518 591 & evap, frugs, agesno, tsoil) … … 748 821 endif 749 822 750 if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice) then 823 if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice 824 : .or.forcing_amma) then 751 825 fcoriolis=0.0 ; ug=0. ; vg=0. 752 826 endif … … 813 887 814 888 teta=temp*(pzero/play)**rkappa 889 ! 890 !--------------------------------------------------------------------- 891 ! Nudge soil temperature if requested 892 !--------------------------------------------------------------------- 893 894 IF (nudge_tsoil) THEN 895 ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) 896 . -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge) 897 ENDIF 815 898 816 899 !---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
