C%%%%%%%%%%%%%%%%%% Abderrahmane 01 06 2005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE RFORCING ( m,rjour,read_oz, e tfi,qfi,wo,wo_dl,wo_ref,wo_dl_ref, e pplayfi,paprsfi,cldfrafi,cldliqfi, e tsolfi,albsfi,rlatfi, s t_seri, s tps_sw_ref,tps_sw_ref0,tps_lw_ref,tps_lw_ref0, s tps_sw_ini,tps_sw_ini0,tps_lw_ini,tps_lw_ini0, s d_tps_sw_ini,d_tps_sw_ini0,d_tps_lw_ini,d_tps_lw_ini0, s d_tps_sw_adj,d_tps_sw_adj0,d_tps_lw_adj,d_tps_lw_adj0, s toa_sw_ref,toa_sw_ref0,toa_lw_ref,toa_lw_ref0, s toa_sw_ini,toa_sw_ini0,toa_lw_ini,toa_lw_ini0, s d_toa_sw_ini,d_toa_sw_ini0,d_toa_lw_ini,d_toa_lw_ini0, s d_toa_sw_adj,d_toa_sw_adj0,d_toa_lw_adj,d_toa_lw_adj0, s srf_sw_ref,srf_sw_ref0,srf_lw_ref,srf_lw_ref0, s d_srf_sw_adj,d_srf_sw_adj0,d_srf_lw_adj,d_srf_lw_adj0, s dHrad_dT,bilq_ref,bilq_ini,bilq_adj, s heat_ref, heat0_ref, cool_ref, cool0_ref, s d_heat_ini, d_heat0_ini, d_cool_ini, d_cool0_ini, s d_heat_adj, d_heat0_adj, d_cool_adj, d_cool0_adj ) C C output C t_seri: adjusted temperature profile ! ! JLD: 2009-07-13 ! Improvment of the convergence of the temperature adjustment. ! Variable dHrad_dT has been introduced. ! ! JLD: 2009-07-12 ! add the flux and the forcing at the surface C use dimphy use IOIPSL ! use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer use radlwsw_m, only: radlwsw IMPLICIT none #include "clesphys.h" #include "nuage.h" #include "thermcell.h" #include "YOMCST.h" LOGICAL debug PARAMETER (debug=.false.) REAL rjour,dtime PARAMETER (dtime=4.*86400.0) INTEGER i, k, j, itap, m, nbs PARAMETER (nbs=30*2) CJLD PARAMETER (nbs=2) REAL co2_ppm_ref, CH4_ppb_ref, N2O_ppb_ref SAVE co2_ppm_ref, CH4_ppb_ref, N2O_ppb_ref REAL CFC11_ppt_ref, CFC12_ppt_ref SAVE CFC11_ppt_ref, CFC12_ppt_ref REAL rco2_ref REAL rch4_ref REAL RN2O_ref Real rcfc11_ref Real rcfc12_ref REAL solaire_ref SAVE rco2_ref,rch4_ref,RN2O_ref, $ rcfc11_ref,rcfc12_ref,solaire_ref REAL co2_ppm_modif,solaire_modif REAL RCO2_modif,RCH4_modif,RN2O_modif, $ RCFC11_modif,RCFC12_modif REAL CH4_ppb_modif,N2O_ppb_modif, $ CFC11_ppt_modif,CFC12_ppt_modif SAVE RCO2_modif,RCH4_modif,RN2O_modif, $ RCFC11_modif,RCFC12_modif,solaire_modif c initialisation des constantes logical ok_newmicro data ok_newmicro/.true./ save ok_newmicro REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon) REAL xflwp(klon), xfiwp(klon) REAL xflwc(klon,klev), xfiwc(klon,klev) REAL sulfate_ref(klon, klev) ! sulfate aerosol mass concentration [ug m-3] REAL sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value) REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3 REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag REAL re(klon, klev) ! cloud droplet effective radius [um] REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell) Real reliq(klon, klev), reice(klon, klev) REAL tfi(klon,klev),qfi(klon,klev),t_seri(klon,klev) REAL paprsfi(klon,klev+1),pplayfi(klon,klev) integer read_climoz,read_oz REAL wo(klon,klev),wo_dl(klon,klev) REAL wo_ref(klon,klev),wo_dl_ref(klon,klev) REAL wof1(klon,klev,1),wof2(klon,klev,2) REAL wof1_ref(klon,klev,1),wof2_ref(klon,klev,2) real dobson_u ! Dobson unit, in kg m-2 parameter (dobson_u=2.1415e-05) real zmasse(klon, klev) REAL albsfi(klon) REAL tsolfi(klon) Cjld integer l,ig0,m integer l,ig0 REAL rlatfi(klon) REAL cldliqfi(klon,klev),cldfrafi(klon,klev) REAL cldtau(klon,klev), cldemi(klon,klev) REAL heat(klon,klev),heat0(klon,klev), $ cool(klon,klev),cool0(klon,klev) REAL heat_ref(klon,klev),heat0_ref(klon,klev), $ cool_ref(klon,klev),cool0_ref(klon,klev) REAL d_heat_ini(klon,klev), d_heat0_ini(klon,klev), $ d_cool_ini(klon,klev), d_cool0_ini(klon,klev) REAL d_heat_adj(klon,klev), d_heat0_adj(klon,klev), $ d_cool_adj(klon,klev), d_cool0_adj(klon,klev) REAL topsw(klon),toplw(klon),solsw(klon),sollw(klon) REAL topsw0(klon),toplw0(klon),solsw0(klon),sollw0(klon) REAL sollwdown(klon),radsol(klon),albpla(klon) REAL flx_lw_up(klon,klev+1),flx_lw_dn(klon,klev+1) REAL flx_lw_up0(klon,klev+1),flx_lw_dn0(klon,klev+1) ! flux LW ciel clair up REAL flx_sw_up(klon,klev+1),flx_sw_dn(klon,klev+1) REAL flx_sw_up0(klon,klev+1),flx_sw_dn0(klon,klev+1) ! flux SW ciel clair up REAL fdh(klon,klev) ! fixed dynamical heating REAL dHrad_dT(klon,klev) ! derivative the radiative heating with temperature c abderrahmane c aerosols real topswad(klon),solswad(klon),topswai(klon),solswai(klon) ! output: aerosol direct forcing at TOA and surface real tau_ae(klon,klev,9,2),piz_ae(klon,klev,9,2) real cg_ae(klon,klev,9,2) ! aerosol optical properties (see aeropt.F) c ! (i.e., with a smaller droplet concentrationand thus larger droplet radii) REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula save bl95_b0, bl95_b1 logical ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not SAVE ok_ade, ok_aie REAL aerindex(klon) ! POLDER aerosol index LOGICAL new_aod save new_aod Real qsat(klon,klev),flwc(klon,klev), fiwc(klon,klev) ! Variable pour iflag_rrtm=1 Real topsw_aero(klon,9),topsw0_aero(klon,9) Real solsw_aero(klon,9), solsw0_aero(klon,9) Real topswad0_aero(klon), solswad0_aero(klon) Real topswcf_aero(klon,3), solswcf_aero(klon,3) logical ok_ade_ref, ok_aie_ref, ok_ade_modif, ok_aie_modif save ok_ade_ref, ok_aie_ref, ok_ade_modif, ok_aie_modif c reference profile REAL tps_sw_ref(klon),tps_lw_ref(klon) REAL tps_sw_ref0(klon),tps_lw_ref0(klon) REAL toa_sw_ref(klon),toa_lw_ref(klon) REAL toa_sw_ref0(klon),toa_lw_ref0(klon) REAL srf_sw_ref(klon),srf_lw_ref(klon) REAL srf_sw_ref0(klon),srf_lw_ref0(klon) c c initial values (with modified PARAMETER, but no strato. adjustment) REAL tps_sw_ini(klon),tps_lw_ini(klon) REAL tps_sw_ini0(klon),tps_lw_ini0(klon) REAL toa_sw_ini(klon),toa_lw_ini(klon) REAL toa_sw_ini0(klon),toa_lw_ini0(klon) c initial anomalies REAL d_tps_sw_ini(klon),d_tps_lw_ini(klon) REAL d_tps_sw_ini0(klon),d_tps_lw_ini0(klon) REAL d_toa_sw_ini(klon),d_toa_lw_ini(klon) REAL d_toa_sw_ini0(klon),d_toa_lw_ini0(klon) c c adjusted values (with modified PARAMETER, and strato. adjustment) REAL tps_sw_adj(klon),tps_lw_adj(klon) REAL tps_sw_adj0(klon),tps_lw_adj0(klon) REAL toa_sw_adj(klon),toa_lw_adj(klon) REAL toa_sw_adj0(klon),toa_lw_adj0(klon) REAL srf_sw_adj(klon),srf_lw_adj(klon) REAL srf_sw_adj0(klon),srf_lw_adj0(klon) c adjusted anomalies REAL d_tps_sw_adj(klon),d_tps_lw_adj(klon) REAL d_tps_sw_adj0(klon),d_tps_lw_adj0(klon) REAL d_toa_sw_adj(klon),d_toa_lw_adj(klon) REAL d_toa_sw_adj0(klon),d_toa_lw_adj0(klon) REAL d_srf_sw_adj(klon),d_srf_lw_adj(klon) REAL d_srf_sw_adj0(klon),d_srf_lw_adj0(klon) c c radiative budget of the stratosphere REAL bilq_ref(klon),bilq_ini(klon),bilq_adj(klon) ! REAL dT_tmp, dT_max(nbs),bilq_max(nbs) c INTEGER kpause(klon) REAL ppause(klon) c REAL dist, rmu0(klon),fract(klon) REAL zlongi real solarlong0 ! si on veut imposer une valeur de la longitude solaire data solarlong0/-999.99/ REAL x_ave, x_std, x_min, x_max REAL undef DATA undef/9999./ logical debut,ok_ozone DATA debut/.true./ save debut if (debut) THEN c print*,'Appel de suphel' CALL suphel !! print*,'Lecture physiq.def' !! !Config Desc = Excentricite !valeur AMIP II R_ecc = 0.016715 call getin('R_ecc', R_ecc) !Config Desc = Equinoxe !valeur AMIP II R_peri = 102.7 call getin('R_peri', R_peri) !Config Desc = Inclinaison !valeur AMIP II R_incl = 23.441 call getin('R_incl', R_incl) !Config Desc = Constante solaire en W/m2 !valeur AMIP II solaire = 1365. call getin('solaire', solaire) ! ! Proprietes des nuages ! rad_froid = 35.0 call getin('rad_froid',rad_froid) ! rad_chau1 = 13.0 call getin('rad_chau1',rad_chau1) ! rad_chau2 = 9.0 call getin('rad_chau2',rad_chau2) ! top_height = 3 call getin('top_height',top_height) ! overlap = 3 call getin('overlap',overlap) ! ! Concentration des gaz a effet de serre ! co2_ppm = 348. call getin('co2_ppm', co2_ppm) ! CH4_ppb = 1650. call getin('CH4_ppb', CH4_ppb) ! N2O_ppb=306. call getin('N2O_ppb', N2O_ppb) ! CFC11_ppt = 280. call getin('CFC11_ppt',CFC11_ppt) ! CFC12_ppt = 484. call getin('CFC12_ppt',CFC12_ppt) ! ! Aerosols ! !Config Desc = Aerosol direct effect not ok_ade = .false. ! call getin('ok_ade', ok_ade) ! !Config Desc = Aerosol indirect effect not ok_aie = .false. ! call getin('ok_aie', ok_aie) ! !Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995) bl95_b0 = 1.7 call getin('bl95_b0', bl95_b0) ! !Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995) bl95_b1 = 0.2 call getin('bl95_b1', bl95_b1) ! !Config Desc = Aerosol ! Test sur new_aod. Ce flag permet de retrouver les resultats de l'AR4 ! il n'est utilisable que lors du couplage avec le SO4 seul new_aod = .false. call getin('new_aod', new_aod) ! !read_climoz read_climoz = 2 call getin('read_climoz', read_climoz) if (read_climoz.ne.read_oz) then print*,'Attention read_climoz de config.def s different de celui correspondant a histmth.nc' endif t_glace_min = 258. call getin('t_glace_min',t_glace_min) t_glace_max = 273.13 call getin('t_glace_max',t_glace_max) ! les valeures lues sont prises comme valeures de références solaire_ref = solaire co2_ppm_ref = co2_ppm CH4_ppb_ref = CH4_ppb N2O_ppb_ref = N2O_ppb CFC11_ppt_ref = CFC11_ppt CFC12_ppt_ref = CFC12_ppt ok_ade_ref = ok_ade ok_aie_ref = ok_aie ! RCO2_ref = co2_ppm_ref * 1.0e-06 * 44.011/28.97 RCH4_ref = CH4_ppb_ref * 1.0E-09 * 16.043/28.97 RN2O_ref = N2O_ppb_ref * 1.0E-09 * 44.013/28.97 RCFC11_ref = CFC11_ppt_ref * 1.0E-12 * 137.3686/28.97 RCFC12_ref = CFC12_ppt_ref * 1.0E-12 * 120.9140/28.97 PRINT*, "Valeurs de rad_chau1, rad_chau2, rad_froid :", . rad_chau1, rad_chau2, rad_froid PRINT* PRINT*, "Lecture des valeures modifiées des paramètres" read(*,*)solaire_modif read(*,*)co2_ppm_modif RCO2_modif = co2_ppm_modif * 1.0e-06 * 44.011/28.97 read(*,*)CH4_ppb_modif RCH4_modif = CH4_ppb_modif * 1.0E-09 * 16.043/28.97 read(*,*)N2O_ppb_modif RN2O_modif = N2O_ppb_modif * 1.0E-09 * 44.013/28.97 read(*,*)CFC11_ppt_modif RCFC11_modif=CFC11_ppt_modif* 1.0E-12 * 137.3686/28.97 read(*,*)CFC12_ppt_modif RCFC12_modif = CFC12_ppt_modif * 1.0E-12 * 120.9140/28.97 read(*,*)ok_ozone ok_ade_modif=.false. ok_aie_modif=.false. ! write(*,*)' ##############################################' write(*,*)' Valeures de références , perturbées et différences' write(*,9000)' Constante solaire =',solaire_ref,solaire_modif $ ,solaire_modif-solaire_ref write(*,9000)' co2_ppm =',co2_ppm_ref,co2_ppm_modif $ ,co2_ppm_modif-co2_ppm_ref write(*,9000)' CH4_ppb =',CH4_ppb_ref,CH4_ppb_modif $ ,CH4_ppb_modif-CH4_ppb_ref write(*,9000)' N2O_ppb =',N2O_ppb_ref,N2O_ppb_modif $ ,N2O_ppb_modif-N2O_ppb_ref write(*,9000)' CFC11_ppt=',CFC11_ppt_ref,CFC11_ppt_modif $ ,CFC11_ppt_modif-CFC11_ppt_ref write(*,9000)' CFC12_ppt=',CFC12_ppt_ref,CFC12_ppt_modif $ ,CFC12_ppt_modif-CFC12_ppt_ref write(*,9001)' ok_ade=',ok_ade_ref, ok_ade_modif write(*,9001)' ok_aie=',ok_aie_ref, ok_aie_modif write(*,*) 9000 FORMAT (1x,A20,3(1pES15.6)) 9001 FORMAT (1x,A20,2(L15)) ! Initialisation rayonnrmrnt RRTM ! call iniradia(klon,klev,paprsfi(1,1:klev+1)) ! debut=.false. endif ! calcul de zmasse if (debug) print*,'rg=',rg forall (k=1: klev) zmasse(:, k)=(paprsfi(:, k)-paprsfi(:, k+1))/rg ! Abderrahmane 26 11 2010 ! A revfaire correctement ! Champ pour RRTM (a lire dans histmth ou recalculer) ! Pour l'instant on les met a 0 qsat=0. flwc=0. fiwc=0. ! initialiser 0 les sorties ! tropopause tps_sw_ref=0. tps_sw_ref0=0. tps_lw_ref=0. tps_lw_ref0=0. tps_sw_ini=0. tps_sw_ini0=0. tps_lw_ini=0. tps_lw_ini0=0. d_tps_sw_ini=0. d_tps_sw_ini0=0. d_tps_lw_ini=0. d_tps_lw_ini0=0. d_tps_sw_adj=0. d_tps_sw_adj0=0. d_tps_lw_adj=0. d_tps_lw_adj0=0. ! toa toa_sw_ref=0. toa_sw_ref0=0. toa_lw_ref=0. toa_lw_ref0=0. toa_sw_ini=0. toa_sw_ini0=0. toa_lw_ini=0. toa_lw_ini0=0. d_toa_sw_ini=0. d_toa_sw_ini0=0. d_toa_lw_ini=0. d_toa_lw_ini0=0. d_toa_sw_adj=0. d_toa_sw_adj0=0. d_toa_lw_adj=0. d_toa_lw_adj0=0. ! surface srf_sw_ref=0. srf_sw_ref0=0. srf_lw_ref=0. srf_lw_ref0=0. d_srf_sw_adj=0. d_srf_sw_adj0=0. d_srf_lw_adj=0. d_srf_lw_adj0=0. ! ! IF (debug) print*,'Appel ozonecm' ! wo(:, :, 1) = ozonecm(rlatfi, paprsfi, rjour) IF (debug) print*,'Appel de tpause' CALL tpause(tfi, pplayfi, kpause, ppause) cCalculer epaisseur optique et emmissivite des nuages IF (debug) print*,'Appel de newmicro' CALL newmicro (paprsfi, pplayfi, ok_newmicro, $ tfi, cldliqfi, cldfrafi, cldtau, cldemi, . pch, pcl, pcm, pct, pctlwp, s xflwp, xfiwp, xflwc, xfiwc, e ok_aie, e sulfate_ref, sulfate_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl, reliq, reice) cpour un jour donne, calculer la longitude vraie de la terre c (par rapport au point vernal-21 mars) dans son orbite solaire c calculer aussi la distance terre-soleil (unite astronomique) IF (debug) print*,'Appel de orbite, rjour=',rjour ! choix entre calcul de la longitude solaire vraie ou valeur fixee a ! solarlong0 ! !!! Atention, modif ad-hoc pour aquaplanette, quand on tourne un seul mois !!!! ! IF ( (m .eq. 1) .and. (rjour. eq.75.) ) solarlong0 = 0. ! IF (debug) print*,'m,rjour,solarlong0:',m,rjour,solarlong0 ! if (solarlong0<-999.) then CALL orbite(rjour,zlongi,dist) if (debug) print*,'dist=',dist ! else ! zlongi=solarlong0 ! longitude solaire vraie ! dist=1. ! distance au soleil / moyenne ! write (*,*) 'Attention, orbite impose, solarlong0=',solarlong0 ! endif cCalculer la duree d'ensoleillement pour un jour et la hauteur c du soleil (cosinus de l'angle zinithal) moyenne sur la journee IF (debug) print*,'Appel angle' CALL angle(zlongi,rlatfi,fract,rmu0) IF (debug) print*,'zlongi:',zlongi,' rlatfi:',rlatfi $ ,' fract:',fract,' rmu0:',rmu0 c DO k = 1, klev DO i=1,klon t_seri(i,k) = tfi(i,k) ENDDO ENDDO c c Appel initial du rayonnement (reference) ! ======================================= c IF (debug) print*,'Appel radlwsw no 1' solaire=solaire_ref RCO2=rco2_ref RCH4=rch4_ref RN2O=RN2O_ref rcfc11=rcfc11_ref rcfc12=rcfc12_ref ok_ade=ok_ade_ref ok_aie=ok_aie_ref IF (debug) print*,'solaire, RCO2, ok_ade, ok_aie', e solaire, RCO2, ok_ade, ok_aie IF (debug) print*,'RCH4, RN2O, rcfc11, rcfc12 ', e RCH4, RN2O, rcfc11, rcfc12 ! ! Attention effet des aerosols dir et indir ! Abderrahmane ! A voir ! tau_ae=0.0 piz_ae=0.0 cg_ae=0.0 !Traitement cas ozone if (read_climoz.eq.2) then wof2_ref(:,:,1)=wo_ref(:,:)/dobson_u/1e3*zmasse(:,:)*rmo3/rmd wof2_ref(:,:,2)=wo_dl_ref(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) if (debug) then print*,'paprsfi(1,:)=',paprsfi(1,:) print*,'pplayfi(1,:)=',pplayfi(1,:) print*,'tsolfi(1)=',tsolfi(1) print*,'albsfi(1)=',albsfi(1) print*,'t_seri(1,:)=',t_seri(1,:) print*,'qfi(1,:)=',qfi(1,:) print*,'wof2_ref(1,:,:)=',wof2_ref(1,:,:) print*,'cldfrafi(1,:)=',cldfrafi(1,:) print*,'cldemi(1,:)=',cldemi(1,:) print*,'cldtau(1,:)=',cldtau(1,:) print*,'ok_ade, ok_aie=',ok_ade, ok_aie endif CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof2_ref, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) else wof1_ref(:,:,1)=wo_ref(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof1_ref, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) endif c DO k = 1, klev !diagnostiquer le FDH DO i = 1,klon fdh(i,k) = -(heat(i,k)-cool(i,k)) ENDDO ENDDO ! ! Flux net DO i = 1, klon ! Ciel avec nuages ! SW tps_sw_ref(i)=flx_sw_dn(i,kpause(i))-flx_sw_up(i,kpause(i)) toa_sw_ref(i)=flx_sw_dn(i,klev+1)-flx_sw_up(i,klev+1) srf_sw_ref(i)=flx_sw_dn(i,1)-flx_sw_up(i,1) ! LW tps_lw_ref(i)=-1.*(flx_lw_dn(i,kpause(i))+flx_lw_up(i,kpause(i))) toa_lw_ref(i)=-1.*(flx_lw_dn(i,klev+1)+flx_lw_up(i,klev+1)) srf_lw_ref(i)=-1.*(flx_lw_dn(i,1)+flx_lw_up(i,1)) ! ! ciel clair ! SW tps_sw_ref0(i)=flx_sw_dn0(i,kpause(i))-flx_sw_up0(i,kpause(i)) toa_sw_ref0(i)=flx_sw_dn0(i,klev+1)-flx_sw_up0(i,klev+1) srf_sw_ref0(i)=flx_sw_dn0(i,1)-flx_sw_up0(i,1) ! LW tps_lw_ref0(i)=-1.*(flx_lw_dn0(i,kpause(i)) $ +flx_lw_up0(i,kpause(i))) toa_lw_ref0(i)=-1.*(flx_lw_dn0(i,klev+1)+flx_lw_up0(i,klev+1)) srf_lw_ref0(i)=-1.*(flx_lw_dn0(i,1)+flx_lw_up0(i,1)) ENDDO ! ! reference radiative inbalance of the stratosphere DO i = 1, klon bilq_ref(i) = (toa_sw_ref(i)-tps_sw_ref(i)) $ +(toa_lw_ref(i)-tps_lw_ref(i)) END DO ! Bilan rad dans l'atmosphere: heat_ref(:,:) = heat(:,:) cool_ref(:,:) = cool(:,:) heat0_ref(:,:) = heat0(:,:) cool0_ref(:,:) = cool0(:,:) ! ! End of computation of the reference case ! ! Calcul de la derivee du bilan en fonction de la temperature ! =========================================================== ! ! On incremente la temp de 1 dans la tropopause DO k = 1, klev DO i = 1, klon IF (k.GE.kpause(i)) t_seri(i,k)=t_seri(i,k)+1. ENDDO ENDDO if (read_climoz.eq.2) then wof2_ref(:,:,1)=wo_ref(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) wof2_ref(:,:,2)=wo_dl_ref(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof2_ref, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) else wof1_ref(:,:,1)=wo_ref(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof1_ref, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) endif ! dHrad_dT(:,:)=-1. DO k = 1, klev !diagnostiquer le d FDH / d T DO i = 1,klon IF (k.GE.kpause(i)) $ dHrad_dT(i,k) = (heat(i,k)-cool(i,k)) + fdh(i,k) ENDDO ENDDO ! On remet la valeure initial du profil de temperature DO k = 1, klev DO i=1,klon t_seri(i,k) = tfi(i,k) ENDDO ENDDO ! Fin du calcul de la derivee du bilan en fonction de la temperature ! ! Calcul des flux radiatifs apres perturbation ! ============================================ ! solaire=solaire_modif RCO2=rco2_modif RCH4=rch4_modif RN2O=RN2O_modif rcfc11=rcfc11_modif rcfc12=rcfc12_modif ok_ade=ok_ade_modif ok_aie=ok_aie_modif if(.NOT.ok_ozone)then wo=wo_ref wo_dl=wo_dl_ref endif ! ! On itere pour ajuster la stratosphere DO 88888 itap = 1, nbs IF (debug) print*,'Appel radlwsw no 2' IF (debug) print*,'solaire, RCO2, ok_ade, ok_aie', e solaire, RCO2, ok_ade, ok_aie if (read_climoz.eq.2) then wof2(:,:,1)=wo(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) wof2(:,:,2)=wo_dl(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof2, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) else wof1(:,:,1)=wo(:,:)/(dobson_u*1e3/zmasse(:,:)/rmo3*rmd) CALL radlwsw(dist, rmu0, fract, e paprsfi, pplayfi,tsolfi,albsfi, e albsfi,t_seri,qfi,wof1, e cldfrafi, cldemi, cldtau, e ok_ade, ok_aie, e tau_ae, piz_ae, cg_ae, e cldtaupi, new_aod, N qsat, flwc, fiwc, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, a sollwdown, s topsw0,toplw0,solsw0,sollw0, s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up, a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up, a topswad, solswad, a topswai, solswai, a topswad0_aero, solswad0_aero, N topsw_aero, topsw0_aero, N solsw_aero, solsw0_aero, N topswcf_aero, solswcf_aero ) endif c Ajouter la tendance des rayonnements et le FDH c dT_max(itap) = 0. DO k = 1, klev DO i = 1, klon IF (k.GE.kpause(i)) THEN dT_tmp = (heat(i,k)-cool(i,k)+fdh(i,k)) $ /(86400./dtime-0.5*dHrad_dT(i,k)) ! print *,i,k, dT_tmp,dHrad_dT(i,k) ! dT_tmp = (heat(i,k)-cool(i,k)+fdh(i,k)) t_seri(i,k) = t_seri(i,k) + dT_tmp dT_max(itap) = max(dT_max(itap),abs(dT_tmp)) ENDIF ENDDO ENDDO c IF ( itap .EQ. 1 ) THEN DO i = 1, klon tps_sw_ini(i) = flx_sw_dn(i,kpause(i)) . - flx_sw_up(i,kpause(i)) toa_sw_ini(i) = flx_sw_dn(i,klev+1) . - flx_sw_up(i,klev+1) tps_lw_ini(i) = -1.*(flx_lw_dn(i,kpause(i)) . + flx_lw_up(i,kpause(i))) toa_lw_ini(i) = -1.*(flx_lw_dn(i,klev+1) . + flx_lw_up(i,klev+1)) c difference relative to the reference values d_tps_sw_ini(i) = tps_sw_ini(i)-tps_sw_ref(i) d_toa_sw_ini(i) = toa_sw_ini(i)-toa_sw_ref(i) d_tps_lw_ini(i) = tps_lw_ini(i)-tps_lw_ref(i) d_toa_lw_ini(i) = toa_lw_ini(i)-toa_lw_ref(i) C tps_sw_ini0(i) = flx_sw_dn0(i,kpause(i)) . - flx_sw_up0(i,kpause(i)) toa_sw_ini0(i) = flx_sw_dn0(i,klev+1) . - flx_sw_up0(i,klev+1) tps_lw_ini0(i) = -1.*(flx_lw_dn0(i,kpause(i)) . + flx_lw_up0(i,kpause(i))) toa_lw_ini0(i) = -1.*(flx_lw_dn0(i,klev+1) . + flx_lw_up0(i,klev+1)) c difference relative to the reference values d_tps_sw_ini0(i) = tps_sw_ini0(i)-tps_sw_ref0(i) d_toa_sw_ini0(i) = toa_sw_ini0(i)-toa_sw_ref0(i) d_tps_lw_ini0(i) = tps_lw_ini0(i)-tps_lw_ref0(i) d_toa_lw_ini0(i) = toa_lw_ini0(i)-toa_lw_ref0(i) END DO ! IF (debug) then CALL cstat(klon,toa_sw_ini,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "toa_sw_ini",x_ave,x_std,x_min,x_max CALL cstat(klon,tps_sw_ini,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "tps_sw_ini",x_ave,x_std,x_min,x_max CALL cstat(klon,toa_sw_ini0,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "toa_sw_ini0",x_ave,x_std,x_min,x_max CALL cstat(klon,tps_sw_ini0,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "tps_sw_ini0",x_ave,x_std,x_min,x_max ! CALL cstat(klon,toa_lw_ini,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "toa_lw_ini",x_ave,x_std,x_min,x_max CALL cstat(klon,tps_lw_ini,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "tps_lw_ini",x_ave,x_std,x_min,x_max CALL cstat(klon,toa_lw_ini0,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "toa_lw_ini0",x_ave,x_std,x_min,x_max CALL cstat(klon,tps_lw_ini0,x_ave,x_std,x_min,x_max,undef) WRITE (*,9002) "tps_lw_ini0",x_ave,x_std,x_min,x_max ! END IF 9002 FORMAT (1x,A12,10(1pE13.6)) ! ! initial radiatiave inbalance of the stratosphere DO i = 1, klon bilq_ini(i) = (toa_sw_ini(i)-tps_sw_ini(i)) $ + (toa_lw_ini(i)-tps_lw_ini(i)) END DO ! ! Variation initiale du bilan rad dans l'atmosphere: d_heat_ini(:,:) = heat(:,:) - heat_ref(:,:) d_cool_ini(:,:) = cool(:,:) - cool_ref(:,:) d_heat0_ini(:,:) = heat0(:,:) - heat0_ref(:,:) d_cool0_ini(:,:) = cool0(:,:) - cool0_ref(:,:) ! END IF ! ( itap .EQ. 1 ) c 88888 CONTINUE DO i = 1, klon ! tps_sw_adj(i) = flx_sw_dn(i,kpause(i)) . - flx_sw_up(i,kpause(i)) toa_sw_adj(i) = flx_sw_dn(i,klev+1) . - flx_sw_up(i,klev+1) srf_sw_adj(i) = flx_sw_dn(i,1) . - flx_sw_up(i,1) tps_lw_adj(i) = -1.*(flx_lw_dn(i,kpause(i)) . + flx_lw_up(i,kpause(i))) toa_lw_adj(i) = -1.*(flx_lw_dn(i,klev+1) . + flx_lw_up(i,klev+1)) srf_lw_adj(i) = -1.*(flx_lw_dn(i,1) . + flx_lw_up(i,1)) c difference relative to the reference values d_tps_sw_adj(i) = tps_sw_adj(i)-tps_sw_ref(i) d_toa_sw_adj(i) = toa_sw_adj(i)-toa_sw_ref(i) d_srf_sw_adj(i) = srf_sw_adj(i)-srf_sw_ref(i) d_tps_lw_adj(i) = tps_lw_adj(i)-tps_lw_ref(i) d_toa_lw_adj(i) = toa_lw_adj(i)-toa_lw_ref(i) d_srf_lw_adj(i) = srf_lw_adj(i)-srf_lw_ref(i) C tps_sw_adj0(i) = flx_sw_dn0(i,kpause(i)) . - flx_sw_up0(i,kpause(i)) toa_sw_adj0(i) = flx_sw_dn0(i,klev+1) . - flx_sw_up0(i,klev+1) srf_sw_adj0(i) = flx_sw_dn0(i,1) . - flx_sw_up0(i,1) tps_lw_adj0(i) = -1.*(flx_lw_dn0(i,kpause(i)) . + flx_lw_up0(i,kpause(i))) toa_lw_adj0(i) = -1.*(flx_lw_dn0(i,klev+1) . + flx_lw_up0(i,klev+1)) srf_lw_adj0(i) = -1.*(flx_lw_dn0(i,1) . + flx_lw_up0(i,1)) c difference relative to the reference values d_tps_sw_adj0(i) = tps_sw_adj0(i)-tps_sw_ref0(i) d_toa_sw_adj0(i) = toa_sw_adj0(i)-toa_sw_ref0(i) d_srf_sw_adj0(i) = srf_sw_adj0(i)-srf_sw_ref0(i) d_tps_lw_adj0(i) = tps_lw_adj0(i)-tps_lw_ref0(i) d_toa_lw_adj0(i) = toa_lw_adj0(i)-toa_lw_ref0(i) d_srf_lw_adj0(i) = srf_lw_adj0(i)-srf_lw_ref0(i) ENDDO c adjusted radiatiave inbalance of the stratosphere DO i = 1, klon bilq_adj(i) = (toa_sw_adj(i)-tps_sw_adj(i)) $ + (toa_lw_adj(i)-tps_lw_adj(i)) END DO ! ! Variation ajustee du bilan rad dans l'atmosphere: d_heat_adj(:,:) = heat(:,:) - heat_ref(:,:) d_cool_adj(:,:) = cool(:,:) - cool_ref(:,:) d_heat0_adj(:,:) = heat0(:,:) - heat0_ref(:,:) d_cool0_adj(:,:) = cool0(:,:) - cool0_ref(:,:) C WRITE (*,9010) 'dT_max',(dT_max(itap),itap=1,nbs) 9010 FORMAT (1x,A8,(10E13.6)) c RETURN END c====================================================================== C SUBROUTINE tpause(t, pls, kpause, ppause) use dimphy IMPLICIT none c REAL plsmin PARAMETER (plsmin=35000.0) REAL seuil PARAMETER (seuil=-2.0e-3) ccc PARAMETER (seuil=2.0e-3) c #include "dimensions.h" REAL t(klon,klev) REAL pls(klon,klev) INTEGER kpause(klon) REAL ppause(klon) c LOGICAL deja(klon) INTEGER kmin(klon) INTEGER i, k REAL deltat, deltaz c DO i = 1, klon deja(i) = .FALSE. ENDDO c DO k = 1, klev DO i = 1, klon IF (.NOT.deja(i) .AND. pls(i,k).LE.plsmin) THEN kmin(i) = k deja(i) = .TRUE. ENDIF ENDDO ENDDO c DO i = 1, klon deja(i) = .FALSE. kpause(i) = kmin(i) ppause(i) = (pls(i,kmin(i))+pls(i,kmin(i)-1))/2.0 ENDDO c DO k = 1, klev DO i = 1, klon IF (.NOT.deja(i) .AND. k.GE.kmin(i)) THEN deltat = t(i,k) - t(i,k-1) deltaz = 287.0/9.81*(t(i,k)+t(i,k-1))/(pls(i,k)+pls(i,k-1)) . *(pls(i,k-1)-pls(i,k)) IF (deltat/deltaz .GT. seuil) THEN C JLD, 2009/11/03: on dit que la limite se situe a la maille du dessous C kpause(i) = k kpause(i) = k-1 CJLD ppause(i) = (pls(i,k)+pls(i,k-1))/2.0 ppause(i) = pls(i,k-1) deja(i) = .TRUE. ENDIF ENDIF ENDDO ENDDO c RETURN END C c====================================================================== C SUBROUTINE lirehist(tapename,iim,jjmp1,klev,ngridmx,mois,m,reado, . tsol,paprs,pplay,t,q,cldfra,cldwcon,ozone,ozone_dl, . rlon, rlat, phis, aire, albs, . rlon_1d, rlat_1d, presnivs) IMPLICIT none c CHARACTER*(*) tapename INTEGER iim, jjmp1, klev, mois, m INTEGER ngridmx REAL rlat(ngridmx), rlon(ngridmx), presnivs(klev) REAL phis(ngridmx), aire(ngridmx), albs(ngridmx) REAL tsol(ngridmx) REAL paprs(ngridmx,klev+1) REAL pplay(ngridmx,klev), t(ngridmx,klev), q(ngridmx,klev) REAL cldfra(ngridmx,klev) REAL cldwcon(ngridmx,klev) ! cloud water content (ice+liquide) integer reado REAL ozone(ngridmx,klev),ozone_dl(ngridmx,klev) REAL cldliq(ngridmx,klev),cldice(ngridmx,klev) REAL psol(ngridmx) C LOGICAL debug PARAMETER (debug=.false.) C PARAMETER (debug=.true.) REAL solsdn(ngridmx), solsup(ngridmx) REAL rlat_1d(jjmp1), rlon_1d(iim) c INTEGER debut4(4), epais4(4), debut3(3), epais3(3) C REAL zx_tmp_2d(iim,jjmp1) REAL zx_tmp_3d(iim,jjmp1,klev) c #include "netcdf.inc" c INTEGER nid, varid, ierr INTEGER i, j, k, l, n c IF (debug) PRINT *,'Entree de lirehist' ierr = NF_OPEN (tapename, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme d ouverture du fichier:"//tapename CALL abort ENDIF c c La taille du fichier en longitude: c IF (debug) PRINT *,'lecture lon' ierr = NF_INQ_DIMID(nid,"lon",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de lon" CALL abort ENDIF ierr = NF_INQ_DIMLEN(nid,varid,i) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir la taille de lon" CALL abort ENDIF IF (i.NE.iim) THEN PRINT*, "Dimension de lon fausse:", i, iim CALL abort ENDIF c c La taille du fichier en latitude: c IF (debug) PRINT *,'lecture lat' ierr = NF_INQ_DIMID(nid,"lat",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de lat" CALL abort ENDIF ierr = NF_INQ_DIMLEN(nid,varid,j) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir la taille de lat" CALL abort ENDIF IF (j.NE.jjmp1) THEN PRINT*, "Dimension lat fausse:", j, jjmp1 CALL abort ENDIF c c La taille du fichier en vertical: c IF (debug) PRINT *,'lecture dimension presniv' ierr = NF_INQ_DIMID(nid,"presnivs",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de z" CALL abort ENDIF ierr = NF_INQ_DIMLEN(nid,varid,k) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir la taille de z" CALL abort ENDIF IF (k.NE.klev) THEN PRINT*, "Dimension z fausse:", k, klev CALL abort ENDIF c c La taille du fichier en temps: c IF (debug) PRINT *,'lecture dimension time' ierr = NF_INQ_DIMID(nid,"time_counter",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de time_counter" CALL abort ENDIF ierr = NF_INQ_DIMLEN(nid,varid,l) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir la taille de time_counter" CALL abort ENDIF IF (m.GT.l) THEN PRINT*, "Dimension time_counter fausse:", l, m CALL abort ENDIF IF (m.GT.mois) THEN PRINT*, "Numero mois sup. dim des mois:", m, mois CALL abort ENDIF c IF (debug) PRINT *,'lecture presniv' ierr = NF_INQ_VARID (nid,"presnivs",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur presnivs" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, varid, presnivs) #else ierr = NF_GET_VAR_REAL(nid, varid, presnivs) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire presnivs" CALL abort ENDIF c ---- IF (debug) PRINT *,'lecture lat' ierr = NF_INQ_VARID (nid,"lat",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur lat" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, varid, rlat_1d) #else ierr = NF_GET_VAR_REAL(nid, varid, rlat_1d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire lat" CALL abort ENDIF DO i=1,iim zx_tmp_2d(i,:)=rlat_1d(:) END DO CALL gr_fi_lu(1,ngridmx,iim,jjmp1,rlat,zx_tmp_2d) c ---- IF (debug) PRINT *,'lecture lon' ierr = NF_INQ_VARID (nid,"lon",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur lon" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, varid, rlon_1d) #else ierr = NF_GET_VAR_REAL(nid, varid, rlon_1d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire lon" CALL abort ENDIF DO j=1,jjmp1 zx_tmp_2d(:,j)=rlon_1d(:) END DO CALL gr_fi_lu(1,ngridmx,iim,jjmp1,rlon,zx_tmp_2d) c c Dans certains fichier hist, phis et aire ont une DIMENSION du temps debut3(1) = 1 debut3(2) = 1 debut3(3) =1 epais3(1) = iim epais3(2) = jjmp1 epais3(3) = 1 C ---- IF (debug) PRINT *,'lecture phis' ierr = NF_INQ_VARID (nid,"phis",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur phis" CALL abort ENDIF #ifdef NC_DOUBLE CJLD ierr = NF_GET_VAR_DOUBLE(nid, varid, phis) ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else CJLD ierr = NF_GET_VAR_REAL(nid, varid, phis) ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire phis" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,phis,zx_tmp_2d) C IF (debug) PRINT *,'lecture aire' ierr = NF_INQ_VARID (nid,"aire",varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur aire" CALL abort ENDIF #ifdef NC_DOUBLE CJLD ierr = NF_GET_VAR_DOUBLE(nid, varid, aire) ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else CJLD ierr = NF_GET_VAR_REAL(nid, varid, aire) ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire aire" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,aire,zx_tmp_2d) c ---- debut3(1) = 1 debut3(2) = 1 debut3(3) = m epais3(1) = iim epais3(2) = jjmp1 epais3(3) = 1 c IF (debug) PRINT *,'lecture tsol' ierr = NF_INQ_VARID(nid,"tsol", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de tsol" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire tsol" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,tsol,zx_tmp_2d) c ---- IF (debug) PRINT *,'lecture psol' ierr = NF_INQ_VARID(nid,"psol", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de psol" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire psol" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,psol,zx_tmp_2d) c ---- IF (debug) PRINT *,'lecture SWdnSFC' ierr = NF_INQ_VARID(nid,"SWdnSFC", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de SWdnSFC" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire SWdnSFC" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,solsdn,zx_tmp_2d) c ---- IF (debug) PRINT *,'lecture SWupSFC' ierr = NF_INQ_VARID(nid,"SWupSFC", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de SWupSFC" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire SWupSFC" CALL abort ENDIF CALL gr_fi_lu(1,ngridmx,iim,jjmp1,solsup,zx_tmp_2d) C C calcul de l'aledo DO n=1,ngridmx IF (solsdn(n).LT.1.0e-3) THEN albs(n)=0.8 ELSE albs(n)=solsup(n)/solsdn(n) ENDIF ENDDO c ---- debut4(1) = 1 debut4(2) = 1 debut4(3) = 1 debut4(4) = m epais4(1) = iim epais4(2) = jjmp1 epais4(3) = klev epais4(4) = 1 c IF (debug) PRINT *,'lecture pres' ierr = NF_INQ_VARID(nid,"pres", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de pres" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire pres" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,pplay,zx_tmp_3d) c DO n = 1, ngridmx paprs(n,1) = psol(n) paprs(n,klev+1) = 0.0 ENDDO DO k = 2, klev DO n = 1, ngridmx paprs(n,k) = (pplay(n,k)+pplay(n,k-1))*0.5 ENDDO ENDDO c ---- IF (debug) PRINT *,'lecture temp' ierr = NF_INQ_VARID(nid,"temp", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de temp" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire temp" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,t,zx_tmp_3d) c ---- IF (debug) PRINT *,'lecture ovap' ierr = NF_INQ_VARID(nid,"ovap", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de ovap" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire ovap" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,q,zx_tmp_3d) cozone ---- IF (debug) PRINT *,'ozone' ! ozone(:,:)=0. ierr = NF_INQ_VARID(nid,"ozone", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur d ozone" CALL abort ENDIF ! ELSE #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire ozone" CALL abort ! ENDIF ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,ozone,zx_tmp_3d) cozone_daylight if (reado.eq.2) then IF (debug) PRINT *,'ozone ozone_daylight' ! ozone(:,:)=0. ierr = NF_INQ_VARID(nid,"ozone_daylight", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur d ozone_dl" CALL abort ENDIF ! ELSE #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire ozone_dl" CALL abort ! ENDIF ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,ozone_dl,zx_tmp_3d) endif c ---- IF (debug) PRINT *,'lecture rneb' ierr = NF_INQ_VARID(nid,"rneb", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de rneb" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire rneb" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldfra,zx_tmp_3d) c ---- IF (debug) PRINT *,'lecture lwcon' ierr = NF_INQ_VARID(nid,"lwcon", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de lwcon" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire lwcon" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldliq,zx_tmp_3d) c ---- c Ice water content IF (debug) PRINT *,'lecture iwcon' ierr = NF_INQ_VARID(nid,"iwcon", varid) IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour avoir identificateur de iwcon" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d) #else ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d) #endif IF (ierr.NE.NF_NOERR) THEN PRINT*, "Probleme pour extraire iwcon" CALL abort ENDIF CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldice,zx_tmp_3d) C cldwcon=cldliq+cldice c ierr = NF_CLOSE(nid) c END C C =================================================== C SUBROUTINE gr_fi_lu(nfield,nlon,iim,jjmp1,fi,lu) IMPLICIT none c c Tranformer une variable de la grille de lecture (hist) a c la grille physique C Routine symetrique de gr_fi_ecrit dans physiq c INTEGER nfield,nlon,iim,jjmp1, jjm REAL fi(nlon,nfield), lu(iim*jjmp1,nfield) c INTEGER i, n, ig c jjm = jjmp1 - 1 DO n = 1, nfield fi(1,n) = lu(1,n) fi(nlon,n) = lu(1+jjm*iim,n) DO ig = 1, nlon - 2 fi(1+ig,n)= lu(iim+ig,n) ENDDO ENDDO RETURN END C --------------------------------------------------------------- SUBROUTINE cstat(ndim,t,ave,std,tmin,tmax,undef) C INTEGER ndim, ntot REAL t(ndim) REAL ave,std, undef REAL*8 tsum,tsum2, tave,tstd REAL tmin,tmax REAL min, max C tsum=0. ntot=0 tmin=1.E31 tmax=-1.E31 DO n=1,ndim C print *,t(n),undef IF (t(n).ne.undef) THEN tsum=tsum+dble(t(n)) tmin=min(tmin,t(n)) tmax=max(tmax,t(n)) ntot=ntot+1 END IF END DO C PRINT *,tsum,ntot tave=tsum/dble(ntot) C tsum2=0. DO n=1,ndim IF (t(n).ne.undef) tsum2=tsum2+dble((t(n)-tave)**2) END DO tstd=sqrt(tsum2/dble(ntot)) C ave=sngl(tave) std=sngl(tstd) C END C C --------------------------------------------------------------- SUBROUTINE cminmax(ndim,t,tmin,tmax,undef) C IMPLICIT none INTEGER ndim, ntot, n REAL t(ndim) REAL tmin,tmax, undef REAL min, max C tmin=1.E31 tmax=-1.E31 DO n=1,ndim C print *,t(n),undef IF (t(n).ne.undef) THEN tmin=min(tmin,t(n)) tmax=max(tmax,t(n)) END IF END DO C END C C --------------------------------------------------------------- REAL FUNCTION coefcorel(ndim,x,y) C INTEGER ndim REAL x(ndim),y(ndim) REAL xave, xstd, yave, ystd, covar, xmin, xmax REAL*8 tsum C CALL cstat(ndim, x, xave, xstd, xmin, xmax, undef) CALL cstat(ndim, y, yave, ystd, xmin, xmax, undef) C tsum=0. DO n=1,ndim tsum=tsum+dble((x(n)-xave)*(y(n)-yave)) END DO covar=sngl(tsum/dble(ndim)) coefcorel=covar/(xstd*ystd) C END c------------------------------------------------------------------- SUBROUTINE cwstat(ndim,t,aire,ave,std,tmin,tmax,undef) ! INTEGER ndim, ntot REAL t(ndim),aire(ndim) REAL ave,std, undef REAL*8 tsum,tsum2, tave,tstd REAL*8 asum,aave REAL tmin,tmax REAL min, max ! tsum=0. asum=0. ntot=0 tmin=1.E31 tmax=-1.E31 DO n=1,ndim ! print *,t(n),undef IF (t(n).ne.undef) THEN tsum=tsum+DBLE(t(n)*aire(n)) asum=asum+DBLE(aire(n)) tmin=min(tmin,t(n)) tmax=max(tmax,t(n)) ntot=ntot+1 END IF END DO ! PRINT *,tsum,ntot tave=tsum/asum aave=asum/dble(ntot) ! tsum2=0. DO n=1,ndim IF (t(n).NE.undef) tsum2=tsum2+DBLE(((t(n)-tave)**2)*aire(n)/aave) END DO tstd=sqrt(tsum2/dble(ntot)) ! ave=sngl(tave) std=sngl(tstd) ! return ! END