Changeset 2160 for LMDZ5/branches/testing/libf/phylmd/1DUTILS.h
- Timestamp:
- Nov 28, 2014, 4:36:29 PM (10 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/1DUTILS.h (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2072,2075-2115,2117-2126,2128-2158
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/1DUTILS.h
r2056 r2160 1 #include "../dyn3d/conf_gcm.F "1 #include "../dyn3d/conf_gcm.F90" 2 2 #include "../dyn3d_common/q_sat.F" 3 3 … … 140 140 CALL getin('ok_flux_surf',ok_flux_surf) 141 141 142 !Config Key = ok_prescr_ust 143 !Config Desc = ustar impose ou non 144 !Config Def = false 145 !Config Help = ustar impose ou non 146 ok_prescr_ust = .false. 147 CALL getin('ok_prescr_ust',ok_prescr_ust) 148 142 149 !Config Key = ok_old_disvert 143 150 !Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) … … 355 362 CHARACTER*(*) fichnom 356 363 !Al1 plev tronque pour .nc mais plev(klev+1):=0 357 real :: plev(klon,klev ),play (klon,klev),phi(klon,klev)364 real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev) 358 365 real :: presnivs(klon,klev) 359 366 real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev) … … 438 445 ! Lecture des champs 439 446 ! 440 plev(1,klev+1)=0.441 CALL get_field("plev",plev,found)442 IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'443 447 CALL get_field("play",play,found) 444 448 IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent' … … 457 461 CALL get_field("omega2",omega2,found) 458 462 IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent' 463 plev(1,klev+1)=0. 464 CALL get_field("plev",plev(:,1:klev),found) 465 IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent' 459 466 460 467 Do iq=1,nqtot … … 822 829 INTEGER np,ierr 823 830 REAL pi,x 824 REAL :: scaleheight=8.825 831 826 832 !----------------------------------------------------------------------- … … 925 931 PRINT *, ap 926 932 927 print*,'scaleheight=',scaleheight 928 DO l = 1, llm 929 dpres(l) = bp(l) - bp(l+1) 930 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 931 write(*, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & 932 log(preff/presnivs(l))*scaleheight & 933 , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ & 934 max(ap(l+1)+bp(l+1)*preff, 1.e-10)) 935 ENDDO 936 933 DO l = 1, llm 934 dpres(l) = bp(l) - bp(l+1) 935 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 936 ENDDO 937 937 938 938 PRINT *,' PRESNIVS ' … … 2108 2108 elseif(4000 < zlay(k) .and. zlay(k) < 9000) then 2109 2109 q_rico(k)=1.8 + (0 - 1.8) / & 2110 & ( 10000 - 4000) * (zlay(k) - 4000)2110 & (9000 - 4000) * (zlay(k) - 4000) 2111 2111 else 2112 2112 q_rico(k)=0.0 … … 2420 2420 end 2421 2421 2422 !===================================================================== 2423 SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof & 2424 & ,th_prof,qv_prof,u_prof,v_prof,o3_prof & 2425 & ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof & 2426 & ,th_mod,qv_mod,u_mod,v_mod,o3_mod & 2427 & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) 2428 2429 implicit none 2430 2431 #include "dimensions.h" 2432 2433 !------------------------------------------------------------------------- 2434 ! Vertical interpolation of Dice forcing data onto model levels 2435 !------------------------------------------------------------------------- 2436 2437 integer nlevmax 2438 parameter (nlevmax=41) 2439 integer nlev_dice,mxcalc,nt_dice 2440 2441 real play(llm), plev_prof(nlev_dice) 2442 real th_prof(nlev_dice),qv_prof(nlev_dice) 2443 real u_prof(nlev_dice),v_prof(nlev_dice) 2444 real o3_prof(nlev_dice) 2445 real ht_prof(nlev_dice),hq_prof(nlev_dice) 2446 real hu_prof(nlev_dice),hv_prof(nlev_dice) 2447 real w_prof(nlev_dice),omega_prof(nlev_dice) 2448 2449 real th_mod(llm),qv_mod(llm) 2450 real u_mod(llm),v_mod(llm), o3_mod(llm) 2451 real ht_mod(llm),hq_mod(llm) 2452 real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm) 2453 2454 integer l,k,k1,k2,kp 2455 real aa,frac,frac1,frac2,fact 2456 2457 do l = 1, llm 2458 2459 if (play(l).ge.plev_prof(nlev_dice)) then 2460 2461 mxcalc=l 2462 k1=0 2463 k2=0 2464 2465 if (play(l).le.plev_prof(1)) then 2466 2467 do k = 1, nlev_dice-1 2468 if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then 2469 k1=k 2470 k2=k+1 2471 endif 2472 enddo 2473 2474 if (k1.eq.0 .or. k2.eq.0) then 2475 write(*,*) 'PB! k1, k2 = ',k1,k2 2476 write(*,*) 'l,play(l) = ',l,play(l)/100 2477 do k = 1, nlev_dice-1 2478 write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100 2479 enddo 2480 endif 2481 2482 frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1)) 2483 th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1)) 2484 qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1)) 2485 u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1)) 2486 v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1)) 2487 o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1)) 2488 ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1)) 2489 hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1)) 2490 hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1)) 2491 hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1)) 2492 w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1)) 2493 omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1)) 2494 2495 else !play>plev_prof(1) 2496 2497 k1=1 2498 k2=2 2499 frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2)) 2500 frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2)) 2501 th_mod(l)= frac1*th_prof(k1) - frac2*th_prof(k2) 2502 qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2) 2503 u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2) 2504 v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2) 2505 o3_mod(l)= frac1*o3_prof(k1) - frac2*o3_prof(k2) 2506 ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2) 2507 hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2) 2508 hu_mod(l)= frac1*hu_prof(k1) - frac2*hu_prof(k2) 2509 hv_mod(l)= frac1*hv_prof(k1) - frac2*hv_prof(k2) 2510 w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2) 2511 omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2) 2512 2513 endif ! play.le.plev_prof(1) 2514 2515 else ! above max altitude of forcing file 2516 2517 !jyg 2518 fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg 2519 fact = max(fact,0.) !jyg 2520 fact = exp(-fact) !jyg 2521 th_mod(l)= th_prof(nlev_dice) !jyg 2522 qv_mod(l)= qv_prof(nlev_dice)*fact !jyg 2523 u_mod(l)= u_prof(nlev_dice)*fact !jyg 2524 v_mod(l)= v_prof(nlev_dice)*fact !jyg 2525 o3_mod(l)= o3_prof(nlev_dice)*fact !jyg 2526 ht_mod(l)= ht_prof(nlev_dice) !jyg 2527 hq_mod(l)= hq_prof(nlev_dice)*fact !jyg 2528 hu_mod(l)= hu_prof(nlev_dice) !jyg 2529 hv_mod(l)= hv_prof(nlev_dice) !jyg 2530 w_mod(l)= 0. !jyg 2531 omega_mod(l)= 0. !jyg 2532 2533 endif ! play 2534 2535 enddo ! l 2536 2537 ! do l = 1,llm 2538 ! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ', 2539 ! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) 2540 ! enddo 2541 2542 return 2543 end 2544 2422 2545 !====================================================================== 2423 2546 SUBROUTINE interp_astex_time(day,day1,annee_ref & … … 2652 2775 2653 2776 !====================================================================== 2777 SUBROUTINE interp_dice_time(day,day1,annee_ref & 2778 & ,year_ini_dice,day_ini_dice,nt_dice,dt_dice & 2779 & ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice & 2780 & ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice & 2781 & ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice & 2782 & ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof & 2783 & ,ustar_prof,psurf_prof,ug_prof,vg_prof & 2784 & ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof) 2785 implicit none 2786 2787 !--------------------------------------------------------------------------------------- 2788 ! Time interpolation of a 2D field to the timestep corresponding to day 2789 ! 2790 ! day: current julian day (e.g. 717538.2) 2791 ! day1: first day of the simulation 2792 ! nt_dice: total nb of data in the forcing (e.g. 145 for Dice) 2793 ! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice) 2794 !--------------------------------------------------------------------------------------- 2795 2796 #include "compar1d.h" 2797 2798 ! inputs: 2799 integer annee_ref 2800 integer nt_dice,nlev_dice 2801 integer year_ini_dice 2802 real day, day1,day_ini_dice,dt_dice 2803 real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice) 2804 real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice) 2805 real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice) 2806 real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice) 2807 real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice) 2808 real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice) 2809 ! outputs: 2810 real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof 2811 real ustar_prof,psurf_prof,ug_prof,vg_prof 2812 real ht_prof(nlev_dice),hq_prof(nlev_dice) 2813 real hu_prof(nlev_dice),hv_prof(nlev_dice) 2814 real w_prof(nlev_dice),omega_prof(nlev_dice) 2815 ! local: 2816 integer it_dice1, it_dice2,k 2817 real timeit,time_dice1,time_dice2,frac 2818 2819 2820 if (forcing_type.eq.7) then 2821 ! Check that initial day of the simulation consistent with Dice period: 2822 print *,'annee_ref=',annee_ref 2823 print *,'day1=',day1 2824 print *,'day_ini_dice=',day_ini_dice 2825 if (annee_ref.ne.1999) then 2826 print*,'Pour Dice, annee_ref doit etre 1999' 2827 print*,'Changer annee_ref dans run.def' 2828 stop 2829 endif 2830 if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then 2831 print*,'Dice a debute le 23 Oct 1999 (jour julien=296)' 2832 print*,'Changer dayref dans run.def',day1,day_ini_dice 2833 stop 2834 endif 2835 if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then 2836 print*,'Dice a fini le 25 Oct 1999 (jour julien=298)' 2837 print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice 2838 stop 2839 endif 2840 2841 endif 2842 2843 ! Determine timestep relative to the 1st day of TOGA-COARE: 2844 ! timeit=(day-day1)*86400. 2845 ! if (annee_ref.eq.1992) then 2846 ! timeit=(day-day_ini_dice)*86400. 2847 ! else 2848 ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992 2849 ! endif 2850 timeit=(day-day_ini_dice)*86400 2851 2852 ! Determine the closest observation times: 2853 it_dice1=INT(timeit/dt_dice)+1 2854 it_dice2=it_dice1 + 1 2855 time_dice1=(it_dice1-1)*dt_dice 2856 time_dice2=(it_dice2-1)*dt_dice 2857 2858 if (it_dice1 .ge. nt_dice) then 2859 write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400. 2860 stop 2861 endif 2862 2863 ! time interpolation: 2864 frac=(time_dice2-timeit)/(time_dice2-time_dice1) 2865 frac=max(frac,0.0) 2866 2867 shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1)) 2868 lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1)) 2869 lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1)) 2870 swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1)) 2871 tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1)) 2872 ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1)) 2873 psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1)) 2874 ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1)) 2875 vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1)) 2876 2877 ! print*, 2878 ! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:', 2879 ! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof 2880 2881 do k=1,nlev_dice 2882 ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1)) 2883 hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1)) 2884 hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1)) 2885 hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1)) 2886 w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1)) 2887 omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1)) 2888 enddo 2889 2890 return 2891 END 2892 2893 !====================================================================== 2654 2894 SUBROUTINE interp_armcu_time(day,day1,annee_ref & 2655 2895 & ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu & … … 2916 3156 return 2917 3157 end 2918 !=====================================================================2919 subroutine read_amma(fich_amma,nlevel,ntime &2920 & ,zz,pp,temp,qv,u,v,dw &2921 & ,dt,dq,sens,flat)2922 2923 !program reading forcings of the AMMA case study2924 2925 2926 implicit none2927 2928 #include "netcdf.inc"2929 2930 integer ntime,nlevel2931 character*80 :: fich_amma2932 real*8 zz(nlevel)2933 2934 real*8 temp(nlevel),pp(nlevel)2935 real*8 qv(nlevel),u(nlevel)2936 real*8 v(nlevel)2937 real*8 dw(nlevel,ntime)2938 real*8 dt(nlevel,ntime)2939 real*8 dq(nlevel,ntime)2940 real*8 flat(ntime),sens(ntime)2941 2942 integer nid, ierr2943 integer nbvar3d2944 parameter(nbvar3d=30)2945 integer var3didin(nbvar3d)2946 2947 ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)2948 if (ierr.NE.NF_NOERR) then2949 write(*,*) 'ERROR: Pb opening forcings nc file '2950 write(*,*) NF_STRERROR(ierr)2951 stop ""2952 endif2953 2954 2955 ierr=NF_INQ_VARID(nid,"zz",var3didin(1))2956 if(ierr/=NF_NOERR) then2957 write(*,*) NF_STRERROR(ierr)2958 stop 'lev'2959 endif2960 2961 2962 ierr=NF_INQ_VARID(nid,"temp",var3didin(2))2963 if(ierr/=NF_NOERR) then2964 write(*,*) NF_STRERROR(ierr)2965 stop 'temp'2966 endif2967 2968 ierr=NF_INQ_VARID(nid,"qv",var3didin(3))2969 if(ierr/=NF_NOERR) then2970 write(*,*) NF_STRERROR(ierr)2971 stop 'qv'2972 endif2973 2974 ierr=NF_INQ_VARID(nid,"u",var3didin(4))2975 if(ierr/=NF_NOERR) then2976 write(*,*) NF_STRERROR(ierr)2977 stop 'u'2978 endif2979 2980 ierr=NF_INQ_VARID(nid,"v",var3didin(5))2981 if(ierr/=NF_NOERR) then2982 write(*,*) NF_STRERROR(ierr)2983 stop 'v'2984 endif2985 2986 ierr=NF_INQ_VARID(nid,"dw",var3didin(6))2987 if(ierr/=NF_NOERR) then2988 write(*,*) NF_STRERROR(ierr)2989 stop 'dw'2990 endif2991 2992 ierr=NF_INQ_VARID(nid,"dt",var3didin(7))2993 if(ierr/=NF_NOERR) then2994 write(*,*) NF_STRERROR(ierr)2995 stop 'dt'2996 endif2997 2998 ierr=NF_INQ_VARID(nid,"dq",var3didin(8))2999 if(ierr/=NF_NOERR) then3000 write(*,*) NF_STRERROR(ierr)3001 stop 'dq'3002 endif3003 3004 ierr=NF_INQ_VARID(nid,"sens",var3didin(9))3005 if(ierr/=NF_NOERR) then3006 write(*,*) NF_STRERROR(ierr)3007 stop 'sens'3008 endif3009 3010 ierr=NF_INQ_VARID(nid,"flat",var3didin(10))3011 if(ierr/=NF_NOERR) then3012 write(*,*) NF_STRERROR(ierr)3013 stop 'flat'3014 endif3015 3016 ierr=NF_INQ_VARID(nid,"pp",var3didin(11))3017 if(ierr/=NF_NOERR) then3018 write(*,*) NF_STRERROR(ierr)3019 stop 'pp'3020 endif3021 3022 !dimensions lecture3023 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)3024 3025 #ifdef NC_DOUBLE3026 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)3027 #else3028 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)3029 #endif3030 if(ierr/=NF_NOERR) then3031 write(*,*) NF_STRERROR(ierr)3032 stop "getvarup"3033 endif3034 ! write(*,*)'lecture z ok',zz3035 3036 #ifdef NC_DOUBLE3037 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)3038 #else3039 ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)3040 #endif3041 if(ierr/=NF_NOERR) then3042 write(*,*) NF_STRERROR(ierr)3043 stop "getvarup"3044 endif3045 ! write(*,*)'lecture th ok',temp3046 3047 #ifdef NC_DOUBLE3048 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)3049 #else3050 ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)3051 #endif3052 if(ierr/=NF_NOERR) then3053 write(*,*) NF_STRERROR(ierr)3054 stop "getvarup"3055 endif3056 ! write(*,*)'lecture qv ok',qv3057 3058 #ifdef NC_DOUBLE3059 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)3060 #else3061 ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)3062 #endif3063 if(ierr/=NF_NOERR) then3064 write(*,*) NF_STRERROR(ierr)3065 stop "getvarup"3066 endif3067 ! write(*,*)'lecture u ok',u3068 3069 #ifdef NC_DOUBLE3070 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)3071 #else3072 ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)3073 #endif3074 if(ierr/=NF_NOERR) then3075 write(*,*) NF_STRERROR(ierr)3076 stop "getvarup"3077 endif3078 ! write(*,*)'lecture v ok',v3079 3080 #ifdef NC_DOUBLE3081 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)3082 #else3083 ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)3084 #endif3085 if(ierr/=NF_NOERR) then3086 write(*,*) NF_STRERROR(ierr)3087 stop "getvarup"3088 endif3089 ! write(*,*)'lecture w ok',dw3090 3091 #ifdef NC_DOUBLE3092 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)3093 #else3094 ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)3095 #endif3096 if(ierr/=NF_NOERR) then3097 write(*,*) NF_STRERROR(ierr)3098 stop "getvarup"3099 endif3100 ! write(*,*)'lecture dt ok',dt3101 3102 #ifdef NC_DOUBLE3103 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)3104 #else3105 ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)3106 #endif3107 if(ierr/=NF_NOERR) then3108 write(*,*) NF_STRERROR(ierr)3109 stop "getvarup"3110 endif3111 ! write(*,*)'lecture dq ok',dq3112 3113 #ifdef NC_DOUBLE3114 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)3115 #else3116 ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)3117 #endif3118 if(ierr/=NF_NOERR) then3119 write(*,*) NF_STRERROR(ierr)3120 stop "getvarup"3121 endif3122 ! write(*,*)'lecture sens ok',sens3123 3124 #ifdef NC_DOUBLE3125 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)3126 #else3127 ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)3128 #endif3129 if(ierr/=NF_NOERR) then3130 write(*,*) NF_STRERROR(ierr)3131 stop "getvarup"3132 endif3133 ! write(*,*)'lecture flat ok',flat3134 3135 #ifdef NC_DOUBLE3136 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)3137 #else3138 ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)3139 #endif3140 if(ierr/=NF_NOERR) then3141 write(*,*) NF_STRERROR(ierr)3142 stop "getvarup"3143 endif3144 ! write(*,*)'lecture pp ok',pp3145 3146 return3147 end subroutine read_amma3148 !======================================================================3149 SUBROUTINE interp_amma_time(day,day1,annee_ref &3150 & ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma &3151 & ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma &3152 & ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)3153 implicit none3154 3155 !---------------------------------------------------------------------------------------3156 ! Time interpolation of a 2D field to the timestep corresponding to day3157 !3158 ! day: current julian day (e.g. 717538.2)3159 ! day1: first day of the simulation3160 ! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)3161 ! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)3162 !---------------------------------------------------------------------------------------3163 3164 #include "compar1d.h"3165 3166 ! inputs:3167 integer annee_ref3168 integer nt_amma,nlev_amma3169 integer year_ini_amma3170 real day, day1,day_ini_amma,dt_amma3171 real vitw_amma(nlev_amma,nt_amma)3172 real ht_amma(nlev_amma,nt_amma)3173 real hq_amma(nlev_amma,nt_amma)3174 real lat_amma(nt_amma)3175 real sens_amma(nt_amma)3176 ! outputs:3177 real vitw_prof(nlev_amma)3178 real ht_prof(nlev_amma)3179 real hq_prof(nlev_amma)3180 real lat_prof,sens_prof3181 ! local:3182 integer it_amma1, it_amma2,k3183 real timeit,time_amma1,time_amma2,frac3184 3185 3186 if (forcing_type.eq.6) then3187 ! Check that initial day of the simulation consistent with AMMA case:3188 if (annee_ref.ne.2006) then3189 print*,'Pour AMMA, annee_ref doit etre 2006'3190 print*,'Changer annee_ref dans run.def'3191 stop3192 endif3193 if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then3194 print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma3195 print*,'Changer dayref dans run.def'3196 stop3197 endif3198 if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then3199 print*,'AMMA a fini le 11 juillet'3200 print*,'Changer dayref ou nday dans run.def'3201 stop3202 endif3203 endif3204 3205 ! Determine timestep relative to the 1st day of AMMA:3206 ! timeit=(day-day1)*86400.3207 ! if (annee_ref.eq.1992) then3208 ! timeit=(day-day_ini_toga)*86400.3209 ! else3210 ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 19923211 ! endif3212 timeit=(day-day_ini_amma)*864003213 3214 ! Determine the closest observation times:3215 ! it_amma1=INT(timeit/dt_amma)+13216 ! it_amma2=it_amma1 + 13217 ! time_amma1=(it_amma1-1)*dt_amma3218 ! time_amma2=(it_amma2-1)*dt_amma3219 3220 it_amma1=INT(timeit/dt_amma)+13221 IF (it_amma1 .EQ. nt_amma) THEN3222 it_amma2=it_amma13223 ELSE3224 it_amma2=it_amma1 + 13225 ENDIF3226 time_amma1=(it_amma1-1)*dt_amma3227 time_amma2=(it_amma2-1)*dt_amma3228 3229 if (it_amma1 .gt. nt_amma) then3230 write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: ' &3231 & ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.3232 stop3233 endif3234 3235 ! time interpolation:3236 frac=(time_amma2-timeit)/(time_amma2-time_amma1)3237 frac=max(frac,0.0)3238 3239 lat_prof = lat_amma(it_amma2) &3240 & -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))3241 sens_prof = sens_amma(it_amma2) &3242 & -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))3243 3244 do k=1,nlev_amma3245 vitw_prof(k) = vitw_amma(k,it_amma2) &3246 & -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))3247 ht_prof(k) = ht_amma(k,it_amma2) &3248 & -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))3249 hq_prof(k) = hq_amma(k,it_amma2) &3250 & -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))3251 enddo3252 3253 return3254 END3255 3158 3256 3159 !===================================================================== … … 3516 3419 return 3517 3420 end subroutine read_fire 3421 !===================================================================== 3422 subroutine read_dice(fich_dice,nlevel,ntime & 3423 & ,zz,pres,th,qv,u,v,o3 & 3424 & ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg & 3425 & ,hadvt,hadvq,hadvu,hadvv,w,omega) 3426 3427 !program reading initial profils and forcings of the Dice case study 3428 3429 3430 implicit none 3431 3432 #include "netcdf.inc" 3433 3434 integer ntime,nlevel 3435 integer l,k 3436 character*80 :: fich_dice 3437 real*8 time(ntime) 3438 real*8 zz(nlevel) 3439 3440 real*8 th(nlevel),pres(nlevel) 3441 real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel) 3442 real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime) 3443 real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime) 3444 real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime) 3445 real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime) 3446 3447 integer nid, ierr 3448 integer nbvar3d 3449 parameter(nbvar3d=30) 3450 integer var3didin(nbvar3d) 3451 3452 ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid) 3453 if (ierr.NE.NF_NOERR) then 3454 write(*,*) 'ERROR: Pb opening forcings nc file ' 3455 write(*,*) NF_STRERROR(ierr) 3456 stop "" 3457 endif 3458 3459 3460 ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 3461 if(ierr/=NF_NOERR) then 3462 write(*,*) NF_STRERROR(ierr) 3463 stop 'height' 3464 endif 3465 3466 ierr=NF_INQ_VARID(nid,"pf",var3didin(11)) 3467 if(ierr/=NF_NOERR) then 3468 write(*,*) NF_STRERROR(ierr) 3469 stop 'pf' 3470 endif 3471 3472 ierr=NF_INQ_VARID(nid,"theta",var3didin(12)) 3473 if(ierr/=NF_NOERR) then 3474 write(*,*) NF_STRERROR(ierr) 3475 stop 'theta' 3476 endif 3477 3478 ierr=NF_INQ_VARID(nid,"qv",var3didin(13)) 3479 if(ierr/=NF_NOERR) then 3480 write(*,*) NF_STRERROR(ierr) 3481 stop 'qv' 3482 endif 3483 3484 ierr=NF_INQ_VARID(nid,"u",var3didin(14)) 3485 if(ierr/=NF_NOERR) then 3486 write(*,*) NF_STRERROR(ierr) 3487 stop 'u' 3488 endif 3489 3490 ierr=NF_INQ_VARID(nid,"v",var3didin(15)) 3491 if(ierr/=NF_NOERR) then 3492 write(*,*) NF_STRERROR(ierr) 3493 stop 'v' 3494 endif 3495 3496 ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16)) 3497 if(ierr/=NF_NOERR) then 3498 write(*,*) NF_STRERROR(ierr) 3499 stop 'o3' 3500 endif 3501 3502 ierr=NF_INQ_VARID(nid,"shf",var3didin(2)) 3503 if(ierr/=NF_NOERR) then 3504 write(*,*) NF_STRERROR(ierr) 3505 stop 'shf' 3506 endif 3507 3508 ierr=NF_INQ_VARID(nid,"lhf",var3didin(3)) 3509 if(ierr/=NF_NOERR) then 3510 write(*,*) NF_STRERROR(ierr) 3511 stop 'lhf' 3512 endif 3513 3514 ierr=NF_INQ_VARID(nid,"lwup",var3didin(4)) 3515 if(ierr/=NF_NOERR) then 3516 write(*,*) NF_STRERROR(ierr) 3517 stop 'lwup' 3518 endif 3519 3520 ierr=NF_INQ_VARID(nid,"swup",var3didin(5)) 3521 if(ierr/=NF_NOERR) then 3522 write(*,*) NF_STRERROR(ierr) 3523 stop 'dqtdx' 3524 endif 3525 3526 ierr=NF_INQ_VARID(nid,"Tg",var3didin(6)) 3527 if(ierr/=NF_NOERR) then 3528 write(*,*) NF_STRERROR(ierr) 3529 stop 'Tg' 3530 endif 3531 3532 ierr=NF_INQ_VARID(nid,"ustar",var3didin(7)) 3533 if(ierr/=NF_NOERR) then 3534 write(*,*) NF_STRERROR(ierr) 3535 stop 'ustar' 3536 endif 3537 3538 ierr=NF_INQ_VARID(nid,"psurf",var3didin(8)) 3539 if(ierr/=NF_NOERR) then 3540 write(*,*) NF_STRERROR(ierr) 3541 stop 'psurf' 3542 endif 3543 3544 ierr=NF_INQ_VARID(nid,"Ug",var3didin(9)) 3545 if(ierr/=NF_NOERR) then 3546 write(*,*) NF_STRERROR(ierr) 3547 stop 'Ug' 3548 endif 3549 3550 ierr=NF_INQ_VARID(nid,"Vg",var3didin(10)) 3551 if(ierr/=NF_NOERR) then 3552 write(*,*) NF_STRERROR(ierr) 3553 stop 'Vg' 3554 endif 3555 3556 ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17)) 3557 if(ierr/=NF_NOERR) then 3558 write(*,*) NF_STRERROR(ierr) 3559 stop 'hadvT' 3560 endif 3561 3562 ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18)) 3563 if(ierr/=NF_NOERR) then 3564 write(*,*) NF_STRERROR(ierr) 3565 stop 'hadvq' 3566 endif 3567 3568 ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19)) 3569 if(ierr/=NF_NOERR) then 3570 write(*,*) NF_STRERROR(ierr) 3571 stop 'hadvu' 3572 endif 3573 3574 ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20)) 3575 if(ierr/=NF_NOERR) then 3576 write(*,*) NF_STRERROR(ierr) 3577 stop 'hadvv' 3578 endif 3579 3580 ierr=NF_INQ_VARID(nid,"w",var3didin(21)) 3581 if(ierr/=NF_NOERR) then 3582 write(*,*) NF_STRERROR(ierr) 3583 stop 'w' 3584 endif 3585 3586 ierr=NF_INQ_VARID(nid,"omega",var3didin(22)) 3587 if(ierr/=NF_NOERR) then 3588 write(*,*) NF_STRERROR(ierr) 3589 stop 'omega' 3590 endif 3591 !dimensions lecture 3592 ! call catchaxis(nid,ntime,nlevel,time,z,ierr) 3593 3594 #ifdef NC_DOUBLE 3595 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz) 3596 #else 3597 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz) 3598 #endif 3599 if(ierr/=NF_NOERR) then 3600 write(*,*) NF_STRERROR(ierr) 3601 stop "getvarup" 3602 endif 3603 ! write(*,*)'lecture zz ok',zz 3604 3605 #ifdef NC_DOUBLE 3606 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres) 3607 #else 3608 ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres) 3609 #endif 3610 if(ierr/=NF_NOERR) then 3611 write(*,*) NF_STRERROR(ierr) 3612 stop "getvarup" 3613 endif 3614 ! write(*,*)'lecture pres ok',pres 3615 3616 #ifdef NC_DOUBLE 3617 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th) 3618 #else 3619 ierr = NF_GET_VAR_REAL(nid,var3didin(12),th) 3620 #endif 3621 if(ierr/=NF_NOERR) then 3622 write(*,*) NF_STRERROR(ierr) 3623 stop "getvarup" 3624 endif 3625 ! write(*,*)'lecture th ok',th 3626 3627 #ifdef NC_DOUBLE 3628 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv) 3629 #else 3630 ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv) 3631 #endif 3632 if(ierr/=NF_NOERR) then 3633 write(*,*) NF_STRERROR(ierr) 3634 stop "getvarup" 3635 endif 3636 ! write(*,*)'lecture qv ok',qv 3637 3638 #ifdef NC_DOUBLE 3639 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u) 3640 #else 3641 ierr = NF_GET_VAR_REAL(nid,var3didin(14),u) 3642 #endif 3643 if(ierr/=NF_NOERR) then 3644 write(*,*) NF_STRERROR(ierr) 3645 stop "getvarup" 3646 endif 3647 ! write(*,*)'lecture u ok',u 3648 3649 #ifdef NC_DOUBLE 3650 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v) 3651 #else 3652 ierr = NF_GET_VAR_REAL(nid,var3didin(15),v) 3653 #endif 3654 if(ierr/=NF_NOERR) then 3655 write(*,*) NF_STRERROR(ierr) 3656 stop "getvarup" 3657 endif 3658 ! write(*,*)'lecture v ok',v 3659 3660 #ifdef NC_DOUBLE 3661 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3) 3662 #else 3663 ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3) 3664 #endif 3665 if(ierr/=NF_NOERR) then 3666 write(*,*) NF_STRERROR(ierr) 3667 stop "getvarup" 3668 endif 3669 ! write(*,*)'lecture o3 ok',o3 3670 3671 #ifdef NC_DOUBLE 3672 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf) 3673 #else 3674 ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf) 3675 #endif 3676 if(ierr/=NF_NOERR) then 3677 write(*,*) NF_STRERROR(ierr) 3678 stop "getvarup" 3679 endif 3680 ! write(*,*)'lecture shf ok',shf 3681 3682 #ifdef NC_DOUBLE 3683 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf) 3684 #else 3685 ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf) 3686 #endif 3687 if(ierr/=NF_NOERR) then 3688 write(*,*) NF_STRERROR(ierr) 3689 stop "getvarup" 3690 endif 3691 ! write(*,*)'lecture lhf ok',lhf 3692 3693 #ifdef NC_DOUBLE 3694 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup) 3695 #else 3696 ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup) 3697 #endif 3698 if(ierr/=NF_NOERR) then 3699 write(*,*) NF_STRERROR(ierr) 3700 stop "getvarup" 3701 endif 3702 ! write(*,*)'lecture lwup ok',lwup 3703 3704 #ifdef NC_DOUBLE 3705 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup) 3706 #else 3707 ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup) 3708 #endif 3709 if(ierr/=NF_NOERR) then 3710 write(*,*) NF_STRERROR(ierr) 3711 stop "getvarup" 3712 endif 3713 ! write(*,*)'lecture swup ok',swup 3714 3715 #ifdef NC_DOUBLE 3716 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg) 3717 #else 3718 ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg) 3719 #endif 3720 if(ierr/=NF_NOERR) then 3721 write(*,*) NF_STRERROR(ierr) 3722 stop "getvarup" 3723 endif 3724 ! write(*,*)'lecture tg ok',tg 3725 3726 #ifdef NC_DOUBLE 3727 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar) 3728 #else 3729 ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar) 3730 #endif 3731 if(ierr/=NF_NOERR) then 3732 write(*,*) NF_STRERROR(ierr) 3733 stop "getvarup" 3734 endif 3735 ! write(*,*)'lecture ustar ok',ustar 3736 3737 #ifdef NC_DOUBLE 3738 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf) 3739 #else 3740 ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf) 3741 #endif 3742 if(ierr/=NF_NOERR) then 3743 write(*,*) NF_STRERROR(ierr) 3744 stop "getvarup" 3745 endif 3746 ! write(*,*)'lecture psurf ok',psurf 3747 3748 #ifdef NC_DOUBLE 3749 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug) 3750 #else 3751 ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug) 3752 #endif 3753 if(ierr/=NF_NOERR) then 3754 write(*,*) NF_STRERROR(ierr) 3755 stop "getvarup" 3756 endif 3757 ! write(*,*)'lecture ug ok',ug 3758 3759 #ifdef NC_DOUBLE 3760 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg) 3761 #else 3762 ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg) 3763 #endif 3764 if(ierr/=NF_NOERR) then 3765 write(*,*) NF_STRERROR(ierr) 3766 stop "getvarup" 3767 endif 3768 ! write(*,*)'lecture vg ok',vg 3769 3770 #ifdef NC_DOUBLE 3771 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt) 3772 #else 3773 ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt) 3774 #endif 3775 if(ierr/=NF_NOERR) then 3776 write(*,*) NF_STRERROR(ierr) 3777 stop "getvarup" 3778 endif 3779 ! write(*,*)'lecture hadvt ok',hadvt 3780 3781 #ifdef NC_DOUBLE 3782 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq) 3783 #else 3784 ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq) 3785 #endif 3786 if(ierr/=NF_NOERR) then 3787 write(*,*) NF_STRERROR(ierr) 3788 stop "getvarup" 3789 endif 3790 ! write(*,*)'lecture hadvq ok',hadvq 3791 3792 #ifdef NC_DOUBLE 3793 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu) 3794 #else 3795 ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu) 3796 #endif 3797 if(ierr/=NF_NOERR) then 3798 write(*,*) NF_STRERROR(ierr) 3799 stop "getvarup" 3800 endif 3801 ! write(*,*)'lecture hadvu ok',hadvu 3802 3803 #ifdef NC_DOUBLE 3804 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv) 3805 #else 3806 ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv) 3807 #endif 3808 if(ierr/=NF_NOERR) then 3809 write(*,*) NF_STRERROR(ierr) 3810 stop "getvarup" 3811 endif 3812 ! write(*,*)'lecture hadvv ok',hadvv 3813 3814 #ifdef NC_DOUBLE 3815 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w) 3816 #else 3817 ierr = NF_GET_VAR_REAL(nid,var3didin(21),w) 3818 #endif 3819 if(ierr/=NF_NOERR) then 3820 write(*,*) NF_STRERROR(ierr) 3821 stop "getvarup" 3822 endif 3823 ! write(*,*)'lecture w ok',w 3824 3825 #ifdef NC_DOUBLE 3826 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega) 3827 #else 3828 ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega) 3829 #endif 3830 if(ierr/=NF_NOERR) then 3831 write(*,*) NF_STRERROR(ierr) 3832 stop "getvarup" 3833 endif 3834 ! write(*,*)'lecture omega ok',omega 3835 3836 return 3837 end subroutine read_dice 3838 !===================================================================== 3839 ! Reads CIRC input files 3840 3841 SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza) 3842 3843 parameter (ncm_1=49180) 3844 #include "YOMCST.h" 3845 3846 real albsfc(ncm_1), albsfc_w(ncm_1) 3847 real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), & 3848 reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ) 3849 real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1) 3850 real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ) 3851 real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ) 3852 real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), & 3853 o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ) 3854 ! za= zenital angle 3855 ! sza= cosinus angle zenital 3856 real wavn(ncm_1), ssf(ncm_1),za,sza 3857 integer nlev 3858 3859 3860 ! Open the files 3861 3862 open (11, file='Tsfc_sza_nlev_case.txt', status='old') 3863 open (12, file='level_input_case.txt', status='old') 3864 open (13, file='layer_input_case.txt', status='old') 3865 open (14, file='aerosol_input_case.txt', status='old') 3866 open (15, file='cloud_input_case.txt', status='old') 3867 open (16, file='sfcalbedo_input_case.txt', status='old') 3868 3869 ! Read scalar information 3870 do iskip=1,5 3871 read (11, *) 3872 enddo 3873 read (11, '(i8)') nlev 3874 read (11, '(f10.2)') tsfc 3875 read (11, '(f10.2)') za 3876 read (11, '(f10.4)') sw_dn_toa 3877 sza=cos(za/180.*RPI) 3878 print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI 3879 close(11) 3880 3881 ! Read level information 3882 read (12, *) 3883 do il=1,nlev 3884 read (12, 302) ilev, z(il), p(il), t(il) 3885 z(il)=z(il)*1000. ! z donne en km 3886 p(il)=p(il)*100. ! p donne en mb 3887 enddo 3888 302 format (i8, f8.3, 2f9.2) 3889 close(12) 3890 3891 ! Read layer information (midpoint values) 3892 do iskip=1,3 3893 read (13, *) 3894 enddo 3895 do il=1,nlev-1 3896 read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), & 3897 n2o(il),co(il),ch4(il),o2(il),ccl4(il), & 3898 f11(il),f12(il) 3899 pm(il)=pm(il)*100. 3900 enddo 3901 303 format (i8, 2f9.2, 10(2x,e13.7)) 3902 close(13) 3903 3904 ! Read aerosol layer information 3905 do iskip=1,3 3906 read (14, *) 3907 enddo 3908 read (14, '(f10.2)') aer_alpha 3909 read (14, *) 3910 read (14, *) 3911 do il=1,nlev-1 3912 read (14, 304) ilev, aer_beta(il), waer(il), gaer(il) 3913 enddo 3914 304 format (i8, f9.5, 2f8.3) 3915 close(14) 3916 3917 ! Read cloud information 3918 do iskip=1,3 3919 read (15, *) 3920 enddo 3921 do il=1,nlev-1 3922 read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il) 3923 lwp(il)=lwp(il)/1000. ! lwp donne en g/kg 3924 iwp(il)=iwp(il)/1000. ! iwp donne en g/kg 3925 reliq(il)=reliq(il)/1000000. ! reliq donne en microns 3926 reice(il)=reice(il)/1000000. ! reice donne en microns 3927 enddo 3928 305 format (i8, f8.3, 4f9.2) 3929 close(15) 3930 3931 ! Read surface albedo (weighted & unweighted) and spectral solar irradiance 3932 do iskip=1,6 3933 read (16, *) 3934 enddo 3935 do icm_1=1,ncm_1 3936 read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1) 3937 enddo 3938 306 format(f10.1, 2f12.5, f14.8) 3939 close(16) 3940 3941 return 3942 end subroutine read_circ 3943 !===================================================================== 3944 ! Reads RTMIP input files 3945 3946 SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3) 3947 3948 #include "YOMCST.h" 3949 3950 real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip) 3951 real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1) 3952 integer nlev 3953 3954 3955 ! Open the files 3956 3957 open (11, file='low_resolution_profile.txt', status='old') 3958 3959 ! Read level information 3960 read (11, *) 3961 do il=1,nlev_rtmip 3962 read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il) 3963 enddo 3964 do il=1,nlev_rtmip 3965 play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb 3966 temp(il)=t(nlev_rtmip-il+1) 3967 ovap(il)=h2o(nlev_rtmip-il+1) 3968 oz(il)=o3(nlev_rtmip-il+1) 3969 enddo 3970 do il=1,39 3971 plev(il)=play(il)+(play(il+1)-play(il))/2. 3972 print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il) 3973 enddo 3974 plev(41)=101300. 3975 302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6) 3976 close(12) 3977 3978 return 3979 end subroutine read_rtmip 3980 !===================================================================== 3981
Note: See TracChangeset
for help on using the changeset viewer.
