1 | subroutine testconserv(ngrid,nlayer,nq,pq,pdq,pqs,pdqs, & |
---|
2 | igcm1,igcm2,ptimestep, & |
---|
3 | pplev,zdq,zdqs,car1,car2) |
---|
4 | |
---|
5 | use comcstfi_mod, only: pi, g |
---|
6 | use geometry_mod, only: latitude, longitude, cell_area |
---|
7 | USE tracer_h, only: igcm_co_ice, igcm_ch4_ice |
---|
8 | USE comgeomfi_h, only: totarea,totarea_planet |
---|
9 | implicit none |
---|
10 | |
---|
11 | !================================================================== |
---|
12 | ! Purpose |
---|
13 | ! ------- |
---|
14 | ! Test conservation of tracers |
---|
15 | ! |
---|
16 | ! Inputs |
---|
17 | ! ------ |
---|
18 | ! ngrid Number of vertical columns |
---|
19 | ! nlayer Number of layers |
---|
20 | ! pplev(ngrid,nlayer+1) Pressure levels |
---|
21 | ! zdq |
---|
22 | ! |
---|
23 | ! Outputs |
---|
24 | ! ------- |
---|
25 | ! |
---|
26 | ! Both |
---|
27 | ! ---- |
---|
28 | ! |
---|
29 | ! Authors |
---|
30 | ! ------- |
---|
31 | ! Tanguy Bertrand 2015 |
---|
32 | ! |
---|
33 | !================================================================== |
---|
34 | |
---|
35 | #include "dimensions.h" |
---|
36 | |
---|
37 | !----------------------------------------------------------------------- |
---|
38 | ! Arguments |
---|
39 | |
---|
40 | INTEGER ngrid, nlayer, nq |
---|
41 | INTEGER igcm1, igcm2 |
---|
42 | REAL pq(ngrid,nlayer,nq) |
---|
43 | REAL pdq(ngrid,nlayer,nq) |
---|
44 | REAL zdq(ngrid,nlayer,nq) |
---|
45 | REAL pqs(ngrid,nq) |
---|
46 | REAL pdqs(ngrid,nq) |
---|
47 | REAL zdqs(ngrid,nq) |
---|
48 | REAL ptimestep |
---|
49 | REAL pplev(ngrid,nlayer+1) |
---|
50 | character(len=3) :: car1 |
---|
51 | character(len=7) :: car2 |
---|
52 | |
---|
53 | ! LOCAL VARIABLES |
---|
54 | INTEGER l,ig,iq |
---|
55 | REAL masse |
---|
56 | REAL pqc(ngrid,nlayer,nq) |
---|
57 | REAL pqcs(ngrid,nq) |
---|
58 | |
---|
59 | ! OUTPUT |
---|
60 | REAL dWtot |
---|
61 | REAL dWtots |
---|
62 | REAL Wtot |
---|
63 | REAL Wtots |
---|
64 | REAL nconsMAX |
---|
65 | INTEGER myig |
---|
66 | REAL vdifcncons(ngrid) |
---|
67 | !----------------------------------------------------------------------- |
---|
68 | |
---|
69 | dWtot=0.0 |
---|
70 | dWtots=0.0 |
---|
71 | Wtot=0.0 |
---|
72 | Wtots=0.0 |
---|
73 | nconsMAX=0.0 |
---|
74 | pqc=pq+pdq*ptimestep |
---|
75 | pqcs=pqs+pdqs*ptimestep |
---|
76 | do ig = 1, ngrid |
---|
77 | vdifcncons(ig)=0.0 |
---|
78 | do l = 1, nlayer |
---|
79 | masse = (pplev(ig,l) - pplev(ig,l+1))/g |
---|
80 | |
---|
81 | iq = igcm1 |
---|
82 | ! sum of atmospheric mass : kg |
---|
83 | Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig) |
---|
84 | dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig) |
---|
85 | ! for each column, total mass lost per sec : kg(tracer) / m2 /s |
---|
86 | vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq) |
---|
87 | |
---|
88 | ! if clouds : |
---|
89 | IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN |
---|
90 | iq = igcm2 |
---|
91 | dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig) |
---|
92 | Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig) |
---|
93 | vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq) |
---|
94 | ENDIF |
---|
95 | |
---|
96 | enddo |
---|
97 | |
---|
98 | iq = igcm1 |
---|
99 | dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig) |
---|
100 | Wtots = Wtots + pqcs(ig,iq)*cell_area(ig) |
---|
101 | vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq) |
---|
102 | |
---|
103 | IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN |
---|
104 | iq = igcm2 |
---|
105 | dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig) |
---|
106 | Wtots = Wtots + pqcs(ig,iq)*cell_area(ig) |
---|
107 | vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq) |
---|
108 | ENDIF |
---|
109 | |
---|
110 | ! vdifcncons is the total amount of material that appear or |
---|
111 | ! disapear per second in the routine |
---|
112 | ! it is the non conservative factor |
---|
113 | |
---|
114 | if(abs(vdifcncons(ig)).gt.abs(nconsMAX))then |
---|
115 | nconsMAX=vdifcncons(ig) |
---|
116 | myig=ig |
---|
117 | endif |
---|
118 | |
---|
119 | enddo |
---|
120 | |
---|
121 | dWtot = dWtot/totarea_planet |
---|
122 | dWtots = dWtots/totarea_planet |
---|
123 | print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2' |
---|
124 | print*,'In ',car2,' surface ',car1,' change =',dWtots,' kg m-2' |
---|
125 | print*,'--> non-cons factor =',dWtot+dWtots,' kg m-2' |
---|
126 | print*,'--> MAX non-cons factor =',nconsMAX,' kg m-2 s-1' |
---|
127 | IF (nconsMAX.gt.0.) then |
---|
128 | print*,'--> obtained at lat/lon=',latitude(myig)*180./pi,longitude(myig)*180./pi |
---|
129 | ENDIF |
---|
130 | print*,'--> Total Mass ',car1,' =',Wtot+Wtots,' kg m-2' |
---|
131 | end subroutine testconserv |
---|
132 | |
---|