source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3d/rforcmain.F @ 1547

Last change on this file since 1547 was 1547, checked in by idelkadi, 13 years ago

Rajout des routines pour le calcul du forcage radiatif "offline" dans LMDZ !
Note technique disponible sur :
​http://www.lmd.jussieu.fr/~idelkadi/Welcome_diagnos.html
  • Property svn:executable set to *
File size: 28.3 KB
Line 
1      PROGRAM rforcmain
2
3      use iophy
4      use ioipsl
5      use dimphy
6      use mod_phys_lmdz_para
7
8      IMPLICIT none
9
10c*******************************************************************************
11c Calcul du forcage radiatif avec ajustement stratospherique
12c A. Idelkadi le 6 mars 2007
13c*******************************************************************************
14C
15! JLD: 2009-07-12
16!    + add the flux and the forcing at the surface
17C JLD 2008-11-29
18C    + propbleme de transfert de la grille de lecture à la grille
19C       physique: debut de re-ecriture
20C    + supression d'une lecture specifique pour l'albedo
21#include "dimensions.h"
22#include "temps.h"
23#include "clesphys.h"
24c======================================================================
25 
26      INTEGER ngridmx
27      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
28 
29c fichier de sortie LMDz contenant les champs qui sera lu pour le calcul
30C des forcages
31      CHARACTER*80 tapename,tapename_ref
32      PARAMETER(tapename_ref="histmth_ref.nc",tapename="histmth.nc")
33      INTEGER nvert ! ID vertical axis
34C
35
36      INTEGER mois
37      PARAMETER (mois=12)
38!       PARAMETER (mois=1) ! changer aussi les DATA de julien
39      REAL dtime ! pas de l'iteration en secondes
40      PARAMETER (dtime=86400.0)
41      real rjour
42      INTEGER jour, julien(mois)
43
44      INTEGER i, k, j, itap, m, n
45
46      REAL rlon_1d(iim),rlat_1d(jjm+1), presl_1d(llm)
47c Variables sur grille physique
48      REAL tfi(ngridmx,llm),qfi(ngridmx,llm),
49     $      tfi_adj(ngridmx,llm)   ! adjusted 3D temp profile
50      REAL paprsfi(ngridmx,llm+1),pplayfi(ngridmx,llm)
51      integer read_climoz ! read ozone climatology
52      Parameter (read_climoz=2)
53      REAL ozone(ngridmx,llm),ozone_ref(ngridmx,llm)
54      Real ozone_daylight(ngridmx,llm),ozone_daylight_ref(ngridmx,llm)
55      REAL rlonfi(ngridmx),rlatfi(ngridmx)
56      REAL cldliqfi(ngridmx,llm),cldfrafi(ngridmx,llm)
57      REAL albsfi(ngridmx)
58      REAL phisfi(ngridmx),airefi(ngridmx),tsolfi(ngridmx)
59
60c     reference profile
61      REAL tps_sw_ref(ngridmx),tps_lw_ref(ngridmx)
62      REAL tps_sw_ref0(ngridmx),tps_lw_ref0(ngridmx)
63      REAL toa_sw_ref(ngridmx),toa_lw_ref(ngridmx)
64      REAL toa_sw_ref0(ngridmx),toa_lw_ref0(ngridmx)
65      REAL srf_sw_ref(ngridmx),srf_lw_ref(ngridmx)
66      REAL srf_sw_ref0(ngridmx),srf_lw_ref0(ngridmx)
67c     initial values (with modofied PARAMETER, but no strato. adjustment)
68      REAL tps_sw_ini(ngridmx),tps_lw_ini(ngridmx)
69      REAL tps_sw_ini0(ngridmx),tps_lw_ini0(ngridmx)
70      REAL toa_sw_ini(ngridmx),toa_lw_ini(ngridmx)
71      REAL toa_sw_ini0(ngridmx),toa_lw_ini0(ngridmx)
72c     initial anomalies
73      REAL d_tps_sw_ini(ngridmx),d_tps_lw_ini(ngridmx)
74      REAL d_tps_sw_ini0(ngridmx),d_tps_lw_ini0(ngridmx)
75      REAL d_toa_sw_ini(ngridmx),d_toa_lw_ini(ngridmx)
76      REAL d_toa_sw_ini0(ngridmx),d_toa_lw_ini0(ngridmx)
77c
78c     adjusted values (with modofied PARAMETER, and strato. adjustment)
79c     adjusted values (with modofied PARAMETER, and strato. adjustment)
80      REAL tps_sw_adj(ngridmx),tps_lw_adj(ngridmx)
81      REAL tps_sw_adj0(ngridmx),tps_lw_adj0(ngridmx)
82      REAL toa_sw_adj(ngridmx),toa_lw_adj(ngridmx)
83      REAL toa_sw_adj0(ngridmx),toa_lw_adj0(ngridmx)
84      REAL srf_sw_adj(ngridmx),srf_lw_adj(ngridmx)
85      REAL srf_sw_adj0(ngridmx),srf_lw_adj0(ngridmx)
86      REAL d_tps_sw_adj(ngridmx),d_tps_lw_adj(ngridmx)
87      REAL d_tps_sw_adj0(ngridmx),d_tps_lw_adj0(ngridmx)
88      REAL d_toa_sw_adj(ngridmx),d_toa_lw_adj(ngridmx)
89      REAL d_toa_sw_adj0(ngridmx),d_toa_lw_adj0(ngridmx)
90      REAL d_srf_sw_adj(ngridmx),d_srf_lw_adj(ngridmx)
91      REAL d_srf_sw_adj0(ngridmx),d_srf_lw_adj0(ngridmx)
92!
93! Bilan radiatif SW (heat) et LW (cool) des differentes couche de l'atmopshre
94      REAL heat_ref(ngridmx,llm), heat0_ref(ngridmx,llm),
95     $     cool_ref(ngridmx,llm), cool0_ref(ngridmx,llm)
96      REAL d_heat_ini(ngridmx,llm), d_heat0_ini(ngridmx,llm),
97     $     d_cool_ini(ngridmx,llm), d_cool0_ini(ngridmx,llm)
98      REAL d_heat_adj(ngridmx,llm), d_heat0_adj(ngridmx,llm),
99     $     d_cool_adj(ngridmx,llm), d_cool0_adj(ngridmx,llm)
100!
101      REAL dHrad_dT(ngridmx,llm) ! derivative the radiative heating with temperature
102      REAL bilq_ref(ngridmx),bilq_ini(ngridmx),bilq_adj(ngridmx)
103
104      DATA (julien(m),m=1,mois) /15, 45, 75, 105, 135, 165,
105     .                           195, 225, 255, 285, 315, 345/
106C      DATA (julien(m),m=1,mois) /75/
107c
108C
109C    Diagnostiques: Valeures moyennes, sur le globe
110C ================================================
111C
112      REAL aire_tot
113C    Valeures mensuelles, moyenne sur le globe
114C
115C     Bilan rad. de la stratosphere
116      REAL bil_str_ref_m(mois), bil_str_ini_m(mois)
117     $    ,bil_str_adj_m(mois)
118!
119C     Flux au sommet de l'atm
120      REAL toa_sw_ref_m(mois), toa_sw_ref0_m(mois)
121     $    , toa_lw_ref_m(mois), toa_lw_ref0_m(mois)
122C     Forcage initial au sommet de l'atm
123      REAL d_toa_sw_ini_m(mois), d_toa_sw_ini0_m(mois)
124     $    ,d_toa_lw_ini_m(mois), d_toa_lw_ini0_m(mois)
125C     Forcage ajuste au sommet de l'atm
126      REAL d_toa_sw_adj_m(mois), d_toa_sw_adj0_m(mois)
127     $    ,d_toa_lw_adj_m(mois), d_toa_lw_adj0_m(mois)
128!
129C     Flux a la tropopause
130      REAL tps_sw_ref_m(mois), tps_sw_ref0_m(mois)
131     $    , tps_lw_ref_m(mois), tps_lw_ref0_m(mois)
132C     Forcage initial a la tropopause
133      REAL d_tps_sw_ini_m(mois), d_tps_sw_ini0_m(mois)
134     $    ,d_tps_lw_ini_m(mois), d_tps_lw_ini0_m(mois)
135C     Forcage ajuste a la tropopause
136      REAL d_tps_sw_adj_m(mois), d_tps_sw_adj0_m(mois)
137     $    ,d_tps_lw_adj_m(mois), d_tps_lw_adj0_m(mois)
138!
139C     Flux a la surface
140      REAL srf_sw_ref_m(mois), srf_sw_ref0_m(mois)
141     $    , srf_lw_ref_m(mois), srf_lw_ref0_m(mois)
142C     Forcage ajuste a la surface
143      REAL d_srf_sw_adj_m(mois), d_srf_sw_adj0_m(mois)
144     $    ,d_srf_lw_adj_m(mois), d_srf_lw_adj0_m(mois)
145!     Bilan radiatif de la stratosphere
146      REAL bilq_ref_m(mois),bilq_ini_m(mois),bilq_adj_m(mois)
147C
148C    Valeures annuelle, moyenne sur le globe
149C ================================================
150C
151C     bilan rad. de la stratosphere
152      REAL bil_str_ref_yr, bil_str_ini_yr, bil_str_adj_yr
153C     Flux au sommet de l'atm
154      REAL toa_sw_ref_yr, toa_sw_ref0_yr
155     $    , toa_lw_ref_yr, toa_lw_ref0_yr
156C     Forcage initial au sommet de l'atm
157      REAL d_toa_sw_ini_yr, d_toa_sw_ini0_yr
158     $    ,d_toa_lw_ini_yr, d_toa_lw_ini0_yr
159C     Forcage ajuste au sommet de l'atm
160      REAL d_toa_sw_adj_yr, d_toa_sw_adj0_yr
161     $    ,d_toa_lw_adj_yr, d_toa_lw_adj0_yr
162!
163C     Flux a la tropopause
164      REAL tps_sw_ref_yr, tps_sw_ref0_yr
165     $    , tps_lw_ref_yr, tps_lw_ref0_yr
166C     Forcage initial a la tropopause
167      REAL d_tps_sw_ini_yr, d_tps_sw_ini0_yr
168     $    ,d_tps_lw_ini_yr, d_tps_lw_ini0_yr
169C     Forcage ajuste a la tropopause
170      REAL d_tps_sw_adj_yr, d_tps_sw_adj0_yr
171     $    ,d_tps_lw_adj_yr, d_tps_lw_adj0_yr
172!
173C     Flux au sommet de l'atm
174      REAL srf_sw_ref_yr, srf_sw_ref0_yr
175     $    , srf_lw_ref_yr, srf_lw_ref0_yr
176C     Forcage ajuste au sommet de l'atm
177      REAL d_srf_sw_adj_yr, d_srf_sw_adj0_yr
178     $    ,d_srf_lw_adj_yr, d_srf_lw_adj0_yr
179C == fin diagnostiques moyenne sur le globe
180C
181c   Initialisation des sorties
182        real zstophy,zstoday,zout,zjulian
183Cjld        integer nhori,nid_day,day_ref
184        integer nhori,nid_day
185        save nid_day
186        INTEGER ndex2d(iim*(jjm+1))
187        REAL zx_tmp_2d(iim,jjm+1)
188        INTEGER ndex3d(iim*(jjm+1)*llm)
189        REAL zx_tmp_3d(iim,jjm+1,llm)
190C
191      LOGICAL aqua_planette
192      PARAMETER (aqua_planette=.false.)
193C
194      LOGICAL debug
195      PARAMETER (debug=.false.)
196C
197        REAL x_ave, x_std, x_min, x_max
198        REAL undef
199        DATA undef/9999./
200c
201        call Init_Phys_lmdz(iim,jjm+1,llm,1,(jjm-1)*iim+2)
202C
203C     Initialisation des variables
204C
205C     Diagnostiques en moyenne sur le globe
206      DO m =1,mois
207C     Flux au sommet de l'atmosphère
208      toa_sw_ref_m(m)  = 0.
209      toa_sw_ref0_m(m) = 0.
210      toa_lw_ref_m(m)  = 0.
211      toa_lw_ref0_m(m) = 0.
212C     Forcage initial au sommet de l'atm
213      d_toa_sw_ini_m(m)  = 0.
214      d_toa_sw_ini0_m(m) = 0.
215      d_toa_lw_ini_m(m)  = 0.
216      d_toa_lw_ini0_m(m) = 0.
217C     Forcage ajuste au sommet de l'atm
218      d_toa_sw_adj_m(m)  = 0.
219      d_toa_sw_adj0_m(m) = 0.
220      d_toa_lw_adj_m(m)  = 0.
221      d_toa_lw_adj0_m(m) = 0.
222C     Flux a la tropopause
223      tps_sw_ref_m(m)  = 0.
224      tps_sw_ref0_m(m) = 0.
225      tps_lw_ref_m(m)  = 0.
226      tps_lw_ref0_m(m) = 0.
227C     Forcage initial a la tropopause
228      d_tps_sw_ini_m(m)  = 0.
229      d_tps_sw_ini0_m(m) = 0.
230      d_tps_lw_ini_m(m)  = 0.
231      d_tps_lw_ini0_m(m) = 0.
232C     Forcage ajuste a la tropopause
233      d_tps_sw_adj_m(m)  = 0.
234      d_tps_sw_adj0_m(m) = 0.
235      d_tps_lw_adj_m(m)  = 0.
236      d_tps_lw_adj0_m(m) = 0.
237C     Flux a la surface
238      srf_sw_ref_m(m)  = 0.
239      srf_sw_ref0_m(m) = 0.
240      srf_lw_ref_m(m)  = 0.
241      srf_lw_ref0_m(m) = 0.
242C     Forcage ajuste a la surface
243      d_srf_sw_adj_m(m)  = 0.
244      d_srf_sw_adj0_m(m) = 0.
245      d_srf_lw_adj_m(m)  = 0.
246      d_srf_lw_adj0_m(m) = 0.
247!     Bilan radiatif de la stratosphere
248      bilq_ref_m(m) = 0.
249      bilq_ini_m(m) = 0.
250      bilq_adj_m(m) = 0.
251C
252      END DO ! boucle sur les mois m
253C     
254C     moyenne annuelle
255C     Flux au sommet de l'atmosphère
256      toa_sw_ref_yr  = 0.
257      toa_sw_ref0_yr = 0.
258      toa_lw_ref_yr  = 0.
259      toa_lw_ref0_yr = 0.
260C     Forcage initial au sommet de l'atm
261      d_toa_sw_ini_yr  = 0.
262      d_toa_sw_ini0_yr = 0.
263      d_toa_lw_ini_yr  = 0.
264      d_toa_lw_ini0_yr = 0.
265C     Forcage ajuste au sommet de l'atm
266      d_toa_sw_adj_yr  = 0.
267      d_toa_sw_adj0_yr = 0.
268      d_toa_lw_adj_yr  = 0.
269      d_toa_lw_adj0_yr = 0.
270C     Flux a la tropopause
271      tps_sw_ref_yr  = 0.
272      tps_sw_ref0_yr = 0.
273      tps_lw_ref_yr  = 0.
274      tps_lw_ref0_yr = 0.
275C     Forcage initial a la tropopause
276      d_tps_sw_ini_yr  = 0.
277      d_tps_sw_ini0_yr = 0.
278      d_tps_lw_ini_yr  = 0.
279      d_tps_lw_ini0_yr = 0.
280C     Forcage ajuste a la tropopause
281      d_tps_sw_adj_yr  = 0.
282      d_tps_sw_adj0_yr = 0.
283      d_tps_lw_adj_yr  = 0.
284      d_tps_lw_adj0_yr = 0.
285C     Flux a la surface
286      srf_sw_ref_yr  = 0.
287      srf_sw_ref0_yr = 0.
288      srf_lw_ref_yr  = 0.
289      srf_lw_ref0_yr = 0.
290C     Forcage ajuste a la surface
291      d_srf_sw_adj_yr  = 0.
292      d_srf_sw_adj0_yr = 0.
293      d_srf_lw_adj_yr  = 0.
294      d_srf_lw_adj0_yr = 0.
295
296c lecture des variables a chaque mois
297C================================================================================
298C*********DEBUT DE LA BOUCLE TEMPORELLE µµµµ************************************µ
299C================================================================================
300        DO 99999 m = 1, mois
301CJLD        DO 99999 m = 1, 2
302
303        WRITE (*,*) 'Mois :',m
304c
305        jour = julien(m)
306        rjour=FLOAT(jour)
307        IF (aqua_planette) rjour=75.
308
309        print*,'Appel de lirehist, rjour=',rjour
310!!! Abderrahmane juin 2011
311! On lit dans le fichier de sortie histmth.nc l'ozone "perturbe"
312       CALL lirehist(tapename_ref,iim,jjm+1,llm,ngridmx,mois,m,
313     .               read_climoz,tsolfi,paprsfi,pplayfi,tfi,qfi,
314     .               cldfrafi,cldliqfi,ozone,ozone_daylight,
315     .               rlonfi, rlatfi, phisfi, airefi, albsfi,
316     .               rlon_1d, rlat_1d, presl_1d)
317
318       CALL lirehist(tapename,iim,jjm+1,llm, ngridmx,mois,m,read_climoz,
319     .       tsolfi,paprsfi,pplayfi,tfi,qfi,cldfrafi,cldliqfi,ozone_ref,
320     .       ozone_daylight_ref,rlonfi, rlatfi, phisfi, airefi, albsfi,
321     .             rlon_1d, rlat_1d, presl_1d)
322       print*,'Apres lirehist'
323C
324       IF (debug)   print*,'Appel ozonecm, rjour=',rjour
325       IF (debug) THEN
326         PRINT *,'klon',klon
327         CALL cstat(klon*klev,ozone,x_ave,x_std,x_min,x_max,undef)
328         WRITE (*,9002) "ozone",x_ave,x_std,x_min,x_max
329       END IF
330! On lit l'ozone directement sur le histmth.nc
331!       CALL ozonecm(rjour, rlatfi, paprsfi, ozone)
332       IF (debug) THEN
333         PRINT *,'klon',klon
334         CALL cstat(klon*klev,ozone,x_ave,x_std,x_min,x_max,undef)
335         WRITE (*,9002) "ozone",x_ave,x_std,x_min,x_max
336       END IF
337
338       IF (debug) THEN
339
340       WRITE (*,*) "ngridmx",ngridmx
341
342       CALL cstat(ngridmx,tsolfi,x_ave,x_std,x_min,x_max,undef)
343       WRITE (*,9002) "tsol",x_ave,x_std,x_min,x_max
344
345       CALL cstat(ngridmx,rlonfi,x_ave,x_std,x_min,x_max,undef)
346       WRITE (*,9002) "rlon",x_ave,x_std,x_min,x_max
347
348       CALL cstat(ngridmx,rlatfi,x_ave,x_std,x_min,x_max,undef)
349       WRITE (*,9002) "rlat",x_ave,x_std,x_min,x_max
350
351       CALL cstat(ngridmx,phisfi,x_ave,x_std,x_min,x_max,undef)
352       WRITE (*,9002) "phis",x_ave,x_std,x_min,x_max
353
354       CALL cstat(ngridmx,airefi,x_ave,x_std,x_min,x_max,undef)
355       WRITE (*,9002) "aire",x_ave,x_std,x_min,x_max
356
357       CALL cstat(ngridmx,albsfi,x_ave,x_std,x_min,x_max,undef)
358       WRITE (*,9002) "albs",x_ave,x_std,x_min,x_max
359
360       CALL cstat(ngridmx*(llm+1),paprsfi,x_ave,x_std,x_min,x_max,undef)
361       WRITE (*,9002) "paprs",x_ave,x_std,x_min,x_max
362
363       CALL cstat(ngridmx*llm,pplayfi,x_ave,x_std,x_min,x_max,undef)
364       WRITE (*,9002) "pplay",x_ave,x_std,x_min,x_max
365
366       CALL cstat(ngridmx*llm,tfi,x_ave,x_std,x_min,x_max,undef)
367       WRITE (*,9002) "temp",x_ave,x_std,x_min,x_max
368
369       CALL cstat(ngridmx*llm,qfi,x_ave,x_std,x_min,x_max,undef)
370       WRITE (*,9002) "ovap",x_ave,x_std,x_min,x_max
371
372       CALL cstat(ngridmx*llm,cldfrafi,x_ave,x_std,x_min,x_max,undef)
373       WRITE (*,9002) "cldfra",x_ave,x_std,x_min,x_max
374
375       CALL cstat(ngridmx*llm,cldliqfi,x_ave,x_std,x_min,x_max,undef)
376       WRITE (*,9002) "cldliq",x_ave,x_std,x_min,x_max
377
378       CALL cstat(ngridmx*llm,ozone,x_ave,x_std,x_min,x_max,undef)
379       WRITE (*,9002) "ozone",x_ave,x_std,x_min,x_max
380
381 9002  FORMAT (1x,A12,10(1pE13.6))
382       END IF
383       if(m.eq.1)then
384#include "ini_histforcing.h"
385        endif
386c
387       heat_ref(:,:) = 0.
388       heat0_ref(:,:) = 0.
389       cool_ref(:,:) = 0.
390       cool0_ref(:,:) = 0.
391!
392       d_heat_ini(:,:) = 0.
393       d_heat0_ini(:,:) = 0.
394       d_cool_ini(:,:) = 0.
395       d_cool0_ini(:,:) = 0.
396!
397       d_heat_adj(:,:) = 0.
398       d_heat0_adj(:,:) = 0.
399       d_cool_adj(:,:) = 0.
400       d_cool0_adj(:,:) = 0.
401!
402       CALL rforcing(m,rjour,read_climoz,
403     $           tfi,qfi,ozone,ozone_daylight,
404     $           ozone_ref,ozone_daylight_ref,
405     $           pplayfi,paprsfi,cldfrafi,cldliqfi,
406     $           tsolfi,albsfi,rlatfi,
407c outputs :
408     $          tfi_adj,
409     $          tps_sw_ref,tps_sw_ref0,tps_lw_ref,tps_lw_ref0,
410     $          tps_sw_ini,tps_sw_ini0,tps_lw_ini,tps_lw_ini0,
411     $  d_tps_sw_ini,d_tps_sw_ini0,d_tps_lw_ini,d_tps_lw_ini0,
412     $  d_tps_sw_adj,d_tps_sw_adj0,d_tps_lw_adj,d_tps_lw_adj0,
413     $          toa_sw_ref,toa_sw_ref0,toa_lw_ref,toa_lw_ref0,
414     $          toa_sw_ini,toa_sw_ini0,toa_lw_ini,toa_lw_ini0,
415     $  d_toa_sw_ini,d_toa_sw_ini0,d_toa_lw_ini,d_toa_lw_ini0,
416     $  d_toa_sw_adj,d_toa_sw_adj0,d_toa_lw_adj,d_toa_lw_adj0,
417     $          srf_sw_ref,srf_sw_ref0,srf_lw_ref,srf_lw_ref0,
418     $  d_srf_sw_adj,d_srf_sw_adj0,d_srf_lw_adj,d_srf_lw_adj0,
419     $  dHrad_dT,bilq_ref,bilq_ini,bilq_adj,
420     $  heat_ref, heat0_ref, cool_ref, cool0_ref,
421     $  d_heat_ini, d_heat0_ini, d_cool_ini, d_cool0_ini,
422     $  d_heat_adj, d_heat0_adj, d_cool_adj, d_cool0_adj )
423C
424      IF (debug) THEN
425        CALL cstat(ngridmx,toa_sw_ref,x_ave,x_std,x_min,x_max,undef)
426        WRITE (*,9002) "toa_sw_ref",x_ave,x_std,x_min,x_max
427        CALL cstat(ngridmx,toa_lw_ref,x_ave,x_std,x_min,x_max,undef)
428        WRITE (*,9002) "toa_lw_ref",x_ave,x_std,x_min,x_max
429        CALL cstat(ngridmx,tps_sw_ref,x_ave,x_std,x_min,x_max,undef)
430        WRITE (*,9002) "tps_sw_ref",x_ave,x_std,x_min,x_max
431        CALL cstat(ngridmx,tps_lw_ref,x_ave,x_std,x_min,x_max,undef)
432        WRITE (*,9002) "tps_lw_ref",x_ave,x_std,x_min,x_max
433        CALL cstat(ngridmx,srf_sw_ref,x_ave,x_std,x_min,x_max,undef)
434        WRITE (*,9002) "srf_sw_ref",x_ave,x_std,x_min,x_max
435        CALL cstat(ngridmx,srf_lw_ref,x_ave,x_std,x_min,x_max,undef)
436        WRITE (*,9002) "srf_lw_ref",x_ave,x_std,x_min,x_max
437!
438        CALL cstat(ngridmx,toa_sw_ref0,x_ave,x_std,x_min,x_max,undef)
439        WRITE (*,9002) "toa_sw_ref0",x_ave,x_std,x_min,x_max
440        CALL cstat(ngridmx,toa_lw_ref0,x_ave,x_std,x_min,x_max,undef)
441        WRITE (*,9002) "toa_lw_ref0",x_ave,x_std,x_min,x_max
442        CALL cwstat(ngridmx,toa_lw_ref0,airefi,x_ave,x_std,x_min,x_max
443     $      ,undef)
444        WRITE (*,9002) "Atoa_lw_ref0",x_ave,x_std,x_min,x_max
445        CALL cstat(ngridmx,tps_sw_ref0,x_ave,x_std,x_min,x_max,undef)
446        WRITE (*,9002) "tps_sw_ref0",x_ave,x_std,x_min,x_max
447        CALL cstat(ngridmx,tps_lw_ref0,x_ave,x_std,x_min,x_max,undef)
448        WRITE (*,9002) "tps_lw_ref0",x_ave,x_std,x_min,x_max
449        CALL cstat(ngridmx,srf_sw_ref0,x_ave,x_std,x_min,x_max,undef)
450        WRITE (*,9002) "srf_sw_ref0",x_ave,x_std,x_min,x_max
451        CALL cstat(ngridmx,srf_lw_ref0,x_ave,x_std,x_min,x_max,undef)
452        WRITE (*,9002) "srf_lw_ref0",x_ave,x_std,x_min,x_max
453!
454        WRITE (*,*)
455!
456        CALL cstat(ngridmx,toa_sw_ini,x_ave,x_std,x_min,x_max,undef)
457        WRITE (*,9002) "toa_sw_ini",x_ave,x_std,x_min,x_max
458        CALL cstat(ngridmx,toa_lw_ini,x_ave,x_std,x_min,x_max,undef)
459        WRITE (*,9002) "toa_lw_ini",x_ave,x_std,x_min,x_max
460        CALL cstat(ngridmx,tps_sw_ini,x_ave,x_std,x_min,x_max,undef)
461        WRITE (*,9002) "tps_sw_ini",x_ave,x_std,x_min,x_max
462        CALL cstat(ngridmx,tps_lw_ini,x_ave,x_std,x_min,x_max,undef)
463        WRITE (*,9002) "tps_lw_ini",x_ave,x_std,x_min,x_max
464!
465        CALL cstat(ngridmx,toa_sw_ini0,x_ave,x_std,x_min,x_max,undef)
466        WRITE (*,9002) "toa_sw_ini0",x_ave,x_std,x_min,x_max
467        CALL cstat(ngridmx,toa_lw_ini0,x_ave,x_std,x_min,x_max,undef)
468        WRITE (*,9002) "toa_lw_ini0",x_ave,x_std,x_min,x_max
469        CALL cwstat(ngridmx,toa_lw_ini0,airefi,x_ave,x_std,x_min,x_max
470     $      ,undef)
471        WRITE (*,9002) "Atoa_lw_ini0",x_ave,x_std,x_min,x_max
472        CALL cstat(ngridmx,tps_sw_ini0,x_ave,x_std,x_min,x_max,undef)
473        WRITE (*,9002) "tps_sw_ini0",x_ave,x_std,x_min,x_max
474        CALL cstat(ngridmx,tps_lw_ini0,x_ave,x_std,x_min,x_max,undef)
475        WRITE (*,9002) "tps_lw_ini0",x_ave,x_std,x_min,x_max
476!
477        CALL cwstat(ngridmx,d_toa_lw_ini0,airefi,x_ave,x_std,x_min,x_max
478     $      ,undef)
479        WRITE (*,9002) "d_toa_lw_ini0",x_ave,x_std,x_min,x_max
480
481      END IF
482
483#include "write_histforcing.h"
484C
485
486C     Diagnostiques en moyenne sur le globe
487      aire_tot = 0.
488      DO n =1, ngridmx
489        aire_tot=aire_tot + airefi(n)
490      END DO
491      DO n =1, ngridmx
492C     Flux au sommet de l'atmosphère
493      toa_sw_ref_m(m)  = toa_sw_ref_m(m)  + toa_sw_ref(n)
494     $      *airefi(n)/aire_tot
495      toa_sw_ref0_m(m) = toa_sw_ref0_m(m) + toa_sw_ref0(n)
496     $      *airefi(n)/aire_tot
497      toa_lw_ref_m(m)  = toa_lw_ref_m(m)  + toa_lw_ref(n)
498     $      *airefi(n)/aire_tot
499      toa_lw_ref0_m(m) = toa_lw_ref0_m(m) + toa_lw_ref0(n)
500     $      *airefi(n)/aire_tot
501C     Forcage initial au sommet de l'atm
502      d_toa_sw_ini_m(m)  = d_toa_sw_ini_m(m)  + d_toa_sw_ini(n)
503     $      *airefi(n)/aire_tot
504      d_toa_sw_ini0_m(m) = d_toa_sw_ini0_m(m) + d_toa_sw_ini0(n)
505     $      *airefi(n)/aire_tot
506      d_toa_lw_ini_m(m)  = d_toa_lw_ini_m(m)  + d_toa_lw_ini(n)
507     $      *airefi(n)/aire_tot
508      d_toa_lw_ini0_m(m) = d_toa_lw_ini0_m(m) + d_toa_lw_ini0(n)
509     $      *airefi(n)/aire_tot
510C     Forcage ajuste au sommet de l'atm
511      d_toa_sw_adj_m(m)  = d_toa_sw_adj_m(m) + d_toa_sw_adj(n)
512     $      *airefi(n)/aire_tot
513      d_toa_sw_adj0_m(m) = d_toa_sw_adj0_m(m) + d_toa_sw_adj0(n)
514     $      *airefi(n)/aire_tot
515      d_toa_lw_adj_m(m)  = d_toa_lw_adj_m(m) + d_toa_lw_adj(n)
516     $      *airefi(n)/aire_tot
517      d_toa_lw_adj0_m(m) = d_toa_lw_adj0_m(m) + d_toa_lw_adj0(n)
518     $      *airefi(n)/aire_tot
519C     Flux a la tropopause
520      tps_sw_ref_m(m)  = tps_sw_ref_m(m)  + tps_sw_ref(n)
521     $      *airefi(n)/aire_tot
522      tps_sw_ref0_m(m) = tps_sw_ref0_m(m) + tps_sw_ref0(n)
523     $      *airefi(n)/aire_tot
524      tps_lw_ref_m(m)  = tps_lw_ref_m(m)  + tps_lw_ref(n)
525     $      *airefi(n)/aire_tot
526      tps_lw_ref0_m(m) = tps_lw_ref0_m(m) + tps_lw_ref0(n)
527     $      *airefi(n)/aire_tot
528C     Forcage initial a la tropopause
529      d_tps_sw_ini_m(m)  = d_tps_sw_ini_m(m)  + d_tps_sw_ini(n)
530     $      *airefi(n)/aire_tot
531      d_tps_sw_ini0_m(m) = d_tps_sw_ini0_m(m) + d_tps_sw_ini0(n)
532     $      *airefi(n)/aire_tot
533      d_tps_lw_ini_m(m)  = d_tps_lw_ini_m(m)  + d_tps_lw_ini(n)
534     $      *airefi(n)/aire_tot
535      d_tps_lw_ini0_m(m) = d_tps_lw_ini0_m(m) + d_tps_lw_ini0(n)
536     $      *airefi(n)/aire_tot
537C     Forcage ajuste a la tropopause
538      d_tps_sw_adj_m(m)  = d_tps_sw_adj_m(m) + d_tps_sw_adj(n)
539     $      *airefi(n)/aire_tot
540      d_tps_sw_adj0_m(m) = d_tps_sw_adj0_m(m) + d_tps_sw_adj0(n)
541     $      *airefi(n)/aire_tot
542      d_tps_lw_adj_m(m)  = d_tps_lw_adj_m(m) + d_tps_lw_adj(n)
543     $      *airefi(n)/aire_tot
544      d_tps_lw_adj0_m(m) = d_tps_lw_adj0_m(m) + d_tps_lw_adj0(n)
545     $      *airefi(n)/aire_tot
546C     Flux a la surface
547      srf_sw_ref_m(m)  = srf_sw_ref_m(m)  + srf_sw_ref(n)
548     $      *airefi(n)/aire_tot
549      srf_sw_ref0_m(m) = srf_sw_ref0_m(m) + srf_sw_ref0(n)
550     $      *airefi(n)/aire_tot
551      srf_lw_ref_m(m)  = srf_lw_ref_m(m)  + srf_lw_ref(n)
552     $      *airefi(n)/aire_tot
553      srf_lw_ref0_m(m) = srf_lw_ref0_m(m) + srf_lw_ref0(n)
554     $      *airefi(n)/aire_tot
555C     Forcage ajuste a la surface
556      d_srf_sw_adj_m(m)  = d_srf_sw_adj_m(m) + d_srf_sw_adj(n)
557     $      *airefi(n)/aire_tot
558      d_srf_sw_adj0_m(m) = d_srf_sw_adj0_m(m) + d_srf_sw_adj0(n)
559     $      *airefi(n)/aire_tot
560      d_srf_lw_adj_m(m)  = d_srf_lw_adj_m(m) + d_srf_lw_adj(n)
561     $      *airefi(n)/aire_tot
562      d_srf_lw_adj0_m(m) = d_srf_lw_adj0_m(m) + d_srf_lw_adj0(n)
563     $      *airefi(n)/aire_tot
564C     Bilan radiatif de la stratosphere
565      bilq_ref_m(m) = bilq_ref_m(m) + bilq_ref(n)
566     $      *airefi(n)/aire_tot
567      bilq_ini_m(m) = bilq_ini_m(m) + bilq_ini(n)
568     $      *airefi(n)/aire_tot
569      bilq_adj_m(m) = bilq_adj_m(m) + bilq_adj(n)
570     $      *airefi(n)/aire_tot
571C
572      END DO
573      CALL cwstat(ngridmx,d_toa_lw_ini0,airefi,x_ave,x_std,x_min,x_max
574     $      ,undef)
575      d_toa_lw_ini0_m(m) = x_ave
576
57799999 CONTINUE
578
579      DO m =1, mois
580C     Flux au sommet de l'atmosphère
581      toa_sw_ref_yr  = toa_sw_ref_yr  + toa_sw_ref_m(m)
582     $      / FLOAT(mois)
583      toa_sw_ref0_yr = toa_sw_ref0_yr + toa_sw_ref0_m(m)
584     $      / FLOAT(mois)
585      toa_lw_ref_yr  = toa_lw_ref_yr  + toa_lw_ref_m(m)
586     $      / FLOAT(mois)
587      toa_lw_ref0_yr = toa_lw_ref0_yr + toa_lw_ref0_m(m)
588     $      / FLOAT(mois)
589C     Forcage initial au sommet de l'atm
590      d_toa_sw_ini_yr  = d_toa_sw_ini_yr  + d_toa_sw_ini_m(m)
591     $      / FLOAT(mois)
592      d_toa_sw_ini0_yr = d_toa_sw_ini0_yr + d_toa_sw_ini0_m(m)
593     $      / FLOAT(mois)
594      d_toa_lw_ini_yr  = d_toa_lw_ini_yr  + d_toa_lw_ini_m(m)
595     $      / FLOAT(mois)
596      d_toa_lw_ini0_yr = d_toa_lw_ini0_yr + d_toa_lw_ini0_m(m)
597     $      / FLOAT(mois)
598C     Forcage ajuste au sommet de l'atm
599      d_toa_sw_adj_yr  = d_toa_sw_adj_yr + d_toa_sw_adj_m(m)
600     $      / FLOAT(mois)
601      d_toa_sw_adj0_yr = d_toa_sw_adj0_yr + d_toa_sw_adj0_m(m)
602     $      / FLOAT(mois)
603      d_toa_lw_adj_yr  = d_toa_lw_adj_yr + d_toa_lw_adj_m(m)
604     $      / FLOAT(mois)
605      d_toa_lw_adj0_yr = d_toa_lw_adj0_yr + d_toa_lw_adj0_m(m)
606     $      / FLOAT(mois)
607C     Flux a la tropopause
608      tps_sw_ref_yr  = tps_sw_ref_yr  + tps_sw_ref_m(m)
609     $      / FLOAT(mois)
610      tps_sw_ref0_yr = tps_sw_ref0_yr + tps_sw_ref0_m(m)
611     $      / FLOAT(mois)
612      tps_lw_ref_yr  = tps_lw_ref_yr  + tps_lw_ref_m(m)
613     $      / FLOAT(mois)
614      tps_lw_ref0_yr = tps_lw_ref0_yr + tps_lw_ref0_m(m)
615     $      / FLOAT(mois)
616C     Forcage initial a la tropopause
617      d_tps_sw_ini_yr  = d_tps_sw_ini_yr  + d_tps_sw_ini_m(m)
618     $      / FLOAT(mois)
619      d_tps_sw_ini0_yr = d_tps_sw_ini0_yr + d_tps_sw_ini0_m(m)
620     $      / FLOAT(mois)
621      d_tps_lw_ini_yr  = d_tps_lw_ini_yr  + d_tps_lw_ini_m(m)
622     $      / FLOAT(mois)
623      d_tps_lw_ini0_yr = d_tps_lw_ini0_yr + d_tps_lw_ini0_m(m)
624     $      / FLOAT(mois)
625C     Forcage ajuste a la tropopause
626      d_tps_sw_adj_yr  = d_tps_sw_adj_yr + d_tps_sw_adj_m(m)
627     $      / FLOAT(mois)
628      d_tps_sw_adj0_yr = d_tps_sw_adj0_yr + d_tps_sw_adj0_m(m)
629     $      / FLOAT(mois)
630      d_tps_lw_adj_yr  = d_tps_lw_adj_yr + d_tps_lw_adj_m(m)
631     $      / FLOAT(mois)
632      d_tps_lw_adj0_yr = d_tps_lw_adj0_yr + d_tps_lw_adj0_m(m)
633     $      / FLOAT(mois)
634C     Flux a la surface
635      srf_sw_ref_yr  = srf_sw_ref_yr  + srf_sw_ref_m(m)
636     $      / FLOAT(mois)
637      srf_sw_ref0_yr = srf_sw_ref0_yr + srf_sw_ref0_m(m)
638     $      / FLOAT(mois)
639      srf_lw_ref_yr  = srf_lw_ref_yr  + srf_lw_ref_m(m)
640     $      / FLOAT(mois)
641      srf_lw_ref0_yr = srf_lw_ref0_yr + srf_lw_ref0_m(m)
642     $      / FLOAT(mois)
643C     Forcage ajuste a la surface
644      d_srf_sw_adj_yr  = d_srf_sw_adj_yr + d_srf_sw_adj_m(m)
645     $      / FLOAT(mois)
646      d_srf_sw_adj0_yr = d_srf_sw_adj0_yr + d_srf_sw_adj0_m(m)
647     $      / FLOAT(mois)
648      d_srf_lw_adj_yr  = d_srf_lw_adj_yr + d_srf_lw_adj_m(m)
649     $      / FLOAT(mois)
650      d_srf_lw_adj0_yr = d_srf_lw_adj0_yr + d_srf_lw_adj0_m(m)
651     $      / FLOAT(mois)
652      END DO
653!
654! Diagnostiques en moyenne mensuelle
655!
656 8001 FORMAT (1x,A10,12(F13.3))
657      WRITE (*,*) 'Stratosphrere radiative budget'
658      WRITE (*,8001) 'Ref :',
659     $     (bilq_ref_m(m),m=1,mois)
660      WRITE (*,8001) 'Ini :',
661     $     (bilq_ini_m(m),m=1,mois)
662      WRITE (*,8001) 'Adj :',
663     $     (bilq_adj_m(m),m=1,mois)
664      WRITE (*,8001) 'Diff :',
665     $     (bilq_adj_m(m)-bilq_ref_m(m),m=1,mois)
666!
667      WRITE (*,*) 'Forcing, adjusted, all sky'
668      WRITE (*,8001) 'TOA SW ',
669     $     (d_toa_sw_adj_m(m),m=1,mois)
670      WRITE (*,8001) 'TOA LW ',
671     $     (d_toa_lw_adj_m(m),m=1,mois)
672      WRITE (*,8001) 'TOA NET',
673     $     ((d_toa_sw_adj_m(m)+d_toa_lw_adj_m(m)),m=1,mois)
674      WRITE (*,8001) 'TPS SW ',
675     $     (d_tps_sw_adj_m(m),m=1,mois)
676      WRITE (*,8001) 'TPS LW ',
677     $     (d_tps_lw_adj_m(m),m=1,mois)
678      WRITE (*,8001) 'TPS NET',
679     $     ((d_tps_sw_adj_m(m)+d_tps_lw_adj_m(m)),m=1,mois)
680      WRITE (*,8001) 'SRF SW ',
681     $     (d_srf_sw_adj_m(m),m=1,mois)
682      WRITE (*,8001) 'SRF LW ',
683     $     (d_srf_lw_adj_m(m),m=1,mois)
684      WRITE (*,8001) 'SRF NET',
685     $     ((d_srf_sw_adj_m(m)+d_srf_lw_adj_m(m)),m=1,mois)
686!
687! Diagnostiques en moyenne annuelle
688!
689 9000 FORMAT (1x,A18,3(A15))
690C 9001 FORMAT (1x,3A9,3(1pEN15.6))
691 9001 FORMAT (1x,A9,A5,A9,3(F13.3))
692
693      WRITE (*,*)
694      WRITE (*,9000) ' ','SW','LW','NET'
695      WRITE (*,9001) 'Bilan','TOA','All sky'
696     $  ,toa_sw_ref_yr, toa_lw_ref_yr, toa_sw_ref_yr + toa_lw_ref_yr
697      WRITE (*,9001) 'Bilan','TPS','All sky'
698     $  ,tps_sw_ref_yr, tps_lw_ref_yr, tps_sw_ref_yr + tps_lw_ref_yr
699      WRITE (*,9001) 'Bilan','SRF','All sky'
700     $  ,srf_sw_ref_yr, srf_lw_ref_yr, srf_sw_ref_yr + srf_lw_ref_yr
701      WRITE (*,*)
702      WRITE (*,9001) 'Ini Forc','TOA','All sky', d_toa_sw_ini_yr
703     $  , d_toa_lw_ini_yr, d_toa_sw_ini_yr + d_toa_lw_ini_yr
704      WRITE (*,9001) 'Adj Forc','TOA','All sky', d_toa_sw_adj_yr
705     $  , d_toa_lw_adj_yr, d_toa_sw_adj_yr + d_toa_lw_adj_yr
706      WRITE (*,*)
707      WRITE (*,9001) 'Ini Forc','TPS','All sky', d_tps_sw_ini_yr
708     $  , d_tps_lw_ini_yr, d_tps_sw_ini_yr + d_tps_lw_ini_yr
709      WRITE (*,9001) 'Adj Forc','TPS','All sky', d_tps_sw_adj_yr
710     $  , d_tps_lw_adj_yr, d_tps_sw_adj_yr + d_tps_lw_adj_yr
711      WRITE (*,*)
712      WRITE (*,9001) 'Adj Forc','SRF','All sky', d_srf_sw_adj_yr
713     $  , d_srf_lw_adj_yr, d_srf_sw_adj_yr + d_srf_lw_adj_yr
714
715      WRITE (*,*) '-----------------------------'
716      WRITE (*,9001) 'Bilan','TOA','Clr sky'
717     $  ,toa_sw_ref0_yr, toa_lw_ref0_yr, toa_sw_ref0_yr + toa_lw_ref0_yr
718      WRITE (*,9001) 'Bilan','TPS','Clr sky'
719     $  ,tps_sw_ref0_yr, tps_lw_ref0_yr, tps_sw_ref0_yr + tps_lw_ref0_yr
720      WRITE (*,9001) 'Bilan','SRF','Clr sky'
721     $  ,srf_sw_ref0_yr, srf_lw_ref0_yr, srf_sw_ref0_yr + srf_lw_ref0_yr
722      WRITE (*,*)
723      WRITE (*,9001) 'Ini Forc','TOA','Clr sky', d_toa_sw_ini0_yr
724     $  , d_toa_lw_ini0_yr, d_toa_sw_ini0_yr + d_toa_lw_ini0_yr
725      WRITE (*,9001) 'Adj Forc','TOA','Clr sky', d_toa_sw_adj0_yr
726     $  , d_toa_lw_adj0_yr, d_toa_sw_adj0_yr + d_toa_lw_adj0_yr
727      WRITE (*,*)
728      WRITE (*,9001) 'Ini Forc','TPS','Clr sky', d_tps_sw_ini0_yr
729     $  , d_tps_lw_ini0_yr, d_tps_sw_ini0_yr + d_tps_lw_ini0_yr
730      WRITE (*,9001) 'Adj Forc','TPS','Clr sky', d_tps_sw_adj0_yr
731     $  , d_tps_lw_adj0_yr, d_tps_sw_adj0_yr + d_tps_lw_adj0_yr
732      WRITE (*,*)
733      WRITE (*,9001) 'Adj Forc','SRF','Clr sky', d_srf_sw_adj0_yr
734     $  , d_srf_lw_adj0_yr, d_srf_sw_adj0_yr + d_srf_lw_adj0_yr
735
736c
737      STOP
738      END
739C
Note: See TracBrowser for help on using the repository browser.