source: trunk/SOURCES/Laure16_files/output_laure16_mod.f90 @ 318

Last change on this file since 318 was 318, checked in by dumas, 4 years ago

Add specific sources for Laure-16 geometry
File size: 15.2 KB
Line 
1module  output_laure16_mod
2
3       USE module3D_phy
4       use bilan_eau_mod
5
6implicit none
7
8real ::  bmean                        !
9real ::  accmean                      ! accumulation moyenne
10real ::  ablmean                      ! ablation moyenne
11real ::  calvmean                     ! moyenne calving
12real ::  ablbordmean                  !
13real ::  ablatotmean
14real ::  bmeltmean                    ! moyenne fusion basale
15real ::  tbmean                       ! temperature basale moyenne
16real ::  tbdotmean                    ! moyenne variation / temps temperature basale
17real ::  vsmean                       ! vitesse de surface moyenne
18!real ::  vsdotmean                    ! moyenne variation / temps vitesse de surface
19real ::  uzsmean   !!!! utilise ?     ! vitesse verticale de surface moyenne
20real ::  uzsdotmean                   ! moyenne variation / temps vitesse verticale de surface
21real ::  uzkmean                      ! moyenne vitesse verticale de surface
22real ::  hdotmean                     ! moyenne derivee / temps de H
23real ::  bdotmean                     ! moyenne bedrock derive / temps
24!real ::  pf0mean                      ! moyenne de PF0
25!real ::  pf1mean                      ! moyenne de PF1
26!real ::  evapmean                     ! moyenne de EVAP
27!real ::  pwmean                       ! moyenne PW
28!real ::  pdfmean                      ! moyenne PDF
29logical,dimension(nx,ny,7) :: mask_cal ! masque regions calotte
30integer, dimension(nx,ny) :: write_mask
31CONTAINS
32!_________________________________________________________________________
33subroutine init_outshort
34
35!ndisp sorite courte tous les ndisp
36NDISP=100
37mask_cal(:,:,:)=.false.
38mask_cal(:,:,1)=.true.
39! calcul d'un masque pour les regions des calottes
40do j=1,ny
41   do i=1,nx
42! -- Laurentide
43      IF( (((xlong(i,j).ge.190).AND.(xlong(i,j).lt.210)).AND. &
44           ((ylat(i,j).ge.50).AND.(ylat(i,j).le.70))).OR.      &
45           (((xlong(i,j).ge.210).AND.(xlong(i,j).lt.220)).AND. &
46           ((ylat(i,j).ge.55).AND.(ylat(i,j).le.75))).OR.      &
47           (((xlong(i,j).ge.220).AND.(xlong(i,j).lt.250)).AND. & 
48           ((ylat(i,j).ge.40).AND.(ylat(i,j).le.85))).OR.      &
49           (((xlong(i,j).ge.250).AND.(xlong(i,j).lt.290)).AND. &
50           ((ylat(i,j).ge.35).AND.(ylat(i,j).le.85))).OR.      &
51           (((xlong(i,j).ge.290).AND.(xlong(i,j).lt.300)).AND. &
52           ((ylat(i,j).ge.35).AND.(ylat(i,j).le.75))).OR.      &
53           (((xlong(i,j).ge.300).AND.(xlong(i,j).le.310)).AND. &
54           ((ylat(i,j).ge.35).AND.(ylat(i,j).lt.60))) ) mask_cal(i,j,2)=.true.
55! -- Groenland
56
57! -- Labrador Sector
58        IF((((xlong(i,j).ge.285).AND.(xlong(i,j).le.300)).AND.  &
59            ((ylat(i,j).ge.35).AND.(ylat(i,j).le.70))).OR.      &
60            (((xlong(i,j).ge.300).AND.(xlong(i,j).lt.310)).AND. &
61            ((ylat(i,j).ge.35).AND.(ylat(i,j).lt.60)))) mask_cal(i,j,4)=.true.
62! -- Keewatin Sector
63        IF(  ((xlong(i,j).gt.240).AND.(xlong(i,j).lt.285)).AND. &
64             ((ylat(i,j).ge.35).AND.(ylat(i,j).le.70)) ) mask_cal(i,j,5)=.true.
65! -- Innuitian Ice Sheet
66        IF( ((xlong(i,j).gt.230).AND.(xlong(i,j).lt.290)).AND. &
67            ((ylat(i,j).gt.70).AND.(ylat(i,j).le.85))) mask_cal(i,j,6)=.true.
68! -- Cordilleran Ice Sheet
69        IF( (((xlong(i,j).ge.190).AND.(xlong(i,j).lt.210)).AND. &
70            ((ylat(i,j).ge.50).AND.(ylat(i,j).le.70))).OR.      &
71            (((xlong(i,j).ge.210).AND.(xlong(i,j).lt.220)).AND. &
72            ((ylat(i,j).ge.55).AND.(ylat(i,j).le.75))).OR.      &
73            (((xlong(i,j).ge.210).AND.(xlong(i,j).lt.220)).AND. &
74             ((ylat(i,j).ge.55).AND.(ylat(i,j).le.75))).OR.     &
75            (((xlong(i,j).ge.220).AND.(xlong(i,j).le.240)).AND. &
76             ((ylat(i,j).ge.40).AND.(ylat(i,j).le.70))) ) mask_cal(i,j,7)=.true.
77! -- Groenland
78        IF( (((xlong(i,j).ge.290).AND.(xlong(i,j).le.310)).AND. &
79            ((ylat(i,j).ge.75).AND.(ylat(i,j).le.85))).OR.      &
80            (((xlong(i,j).ge.300).AND.(xlong(i,j).lt.310)).AND. &
81            ((ylat(i,j).ge.60).AND.(ylat(i,j).le.75))).OR.      &
82!           (((xlong(i,j).ge.310).AND.(xlong(i,j).le.350)).AND. &
83!           ((ylat(i,j).ge.40).AND.(ylat(i,j).le.85))) ) THEN
84            (((xlong(i,j).ge.310).AND.(xlong(i,j).le.345)).AND. &
85            ((ylat(i,j).ge.54).AND.(ylat(i,j).le.85))) ) mask_cal(i,j,3)=.true.
86       enddo
87    enddo
88!        open (UNIT=7777,file='mask_hemin40-02.dat')
89!   do j=1,ny
90!        do i=1,nx
91!        if (mask_cal(i,j,2)) then
92!                write_mask(i,j)=2
93!        else
94!                write_mask(i,j)=0
95!        endif   
96!        write(7777,*) i,j,write_mask(i,j)
97!        enddo     
98!   enddo
99!   close(UNIT=7777)
100!   stop
101
102! ecriture entete avec les colonnes du fichier short :
103WRITE(num_ritz, '(a91)')"! hemin40    : time,isvol,deltavol,inp,isvolf,inf,isacc,isabl,ISABLBORD,ABLATOT,ISCALV,ISBM"
104WRITE(num_ritz, '(a90)')"! colonne    :  1    2      3       4    5     6    7     8      9        10      11    12"
105WRITE(num_ritz,'(a157)')"! calotte    : Laurentide Groenland LabradorSector KeewatinSector InnuitianIS CordilleranIS"
106WRITE(num_ritz,'(a153)')"! isvol      :    13        26          39            52            65"
107WRITE(num_ritz,'(a153)')"! inp        :    14        27          40            53            66"
108WRITE(num_ritz,'(a153)')"! isvolf     :    15        28          41            54            67"
109WRITE(num_ritz,'(a153)')"! inf        :    16        29          42            55            68"
110WRITE(num_ritz,'(a153)')"! isacc      :    17        30          43            56            69"
111WRITE(num_ritz,'(a153)')"! isabl      :    18        31          44            57            70"
112WRITE(num_ritz,'(a153)')"! isablbord  :    19        32          45            58            71"
113WRITE(num_ritz,'(a153)')"! ablatot    :    20        33          46            59            72"
114WRITE(num_ritz,'(a153)')"! iscalv     :    21        34          47            60            73"
115WRITE(num_ritz,'(a153)')"! isbm       :    22        35          48            61            74"
116WRITE(num_ritz,'(a153)')"! itjjamean_ :    23        36          49            62            75"
117WRITE(num_ritz,'(a153)')"! hmean_     :    24        37          50            63            76"
118WRITE(num_ritz,'(a153)')"! Hmax_      :    25        38          51            64            77"
119WRITE(num_ritz,*)
120end subroutine init_outshort
121
122!_________________________________________________________________________
123subroutine shortoutput
124
125! 1_initialization
126!------------------
127      integer KK
128!     integer inpl, INPG, INPF
129!     integer inplab,inpkeew,inpinn,inpcord
130!     integer inpscaN,inpbar,inpkara,inpsib,inpbrit
131     integer inp(7) ! Surface posee (nb de noeuds)
132     integer inf(7) ! surface flottante
133
134!Variables pour sommer     
135      real isvol(7),isvolf(7) ! volume posé et flottants
136      REAL ISCALV(7),ISACC(7),ISBM(7),ISABL(7)
137      REAL ISABLBORD(7),ABLATOT(7),TACC(7),TBM(7)
138      REAL ITJJA(7)
139!      REAL ABLAMEAN
140!moyennes utilisées en output
141      REAL HMAX_(7) , HMEAN_(7) 
142      REAL BMEAN_(7) , ACCMEAN_(7) , ABLMEAN_(7) , ABLBORDMEAN_(7)
143      REAL ABLATOTMEAN_(7) , CALVMEAN_(7) , ITJJAMEAN_(7)
144
145
146
147!     REAL ITJJAMEAN_L, ITJJAMEAN_G, ITJJAMEAN_F ...
148!     REAL ITJJAMEAN_LAB,ITJJAMEAN_KEEW,ITJJAMEAN_INN,ITJJAMEAN_CORD
149!     REAL ITJJAMEAN_SCAN,ITJJAMEAN_BAR,ITJJAMEAN_KARA,ITJJAMEAN_SIB,ITJJAMEAN_BRIT
150
151      REAL HDOTMEAN_G
152!      REAL ABLA(NX,NY)
153      REAL DELTAVOL
154
155
156!     open(unit=4145,file='reg_output_nord.dat')
157
158      BMELTMEAN=0.
159      TBMEAN=0.
160      TBDOTMEAN=0.
161      VSMEAN=0.
162      UZSMEAN=0.
163      UZSDOTMEAN=0.
164      UZKMEAN=0.
165      HDOTMEAN=0.
166      HDOTMEAN_G = 0.
167      BDOTMEAN=0.
168
169      DO kk = 1,13
170        INP(kk) = 0
171        INF(kk) = 0
172        ISVOL(kk) = 0.
173        ISVOLF(kk) = 0.
174        ISBM(kk) = 0.
175        ISACC(kk) = 0.
176        ISCALV(kk) = 0.
177        ISABL(kk) = 0.
178        ISABLBORD(kk) = 0.
179        ABLATOT(kk)= 0.
180        TACC(kk) = 0.
181        TBM(kk) = 0.
182        ITJJA(kk) = 0.
183       
184! nouveau tof mai 2009
185!        where (mask_cal(:,:,kk).and.H(:,:).gt.2..and.flot(:,:)) ISVOLF(kk) = ISVOLF(kk) + H(I,J)
186        isvolf(kk) = sum(H(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and.flot(:,:)))
187        INF(kk) = count(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and.flot(:,:))
188
189        INP(kk) = count(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:))
190        isvol(kk) = sum(H(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:)))
191        isacc(kk) = sum(Acc(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:)))
192        isabl(kk) = sum(Abl(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:)))
193        itjja(kk) = sum(Tjuly(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:)))
194        Hmax_(kk) = maxval(H(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.).and..not.flot(:,:)))
195        Tacc(kk) = sum(Acc(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.)))
196        isablbord(kk) = sum(ablbord(:,:),mask=(mask_cal(:,:,kk)))
197        iscalv(kk) =  sum(calv(:,:),mask=(mask_cal(:,:,kk).and.(H(:,:).gt.2.)))/dt
198
199        ablatot(kk) = isabl(kk) + isablbord(kk)
200        isbm(kk) = isacc(kk)+isabl(kk)+iscalv(kk)+isablbord(kk)
201        Tbm(kk) = Tacc(kk)+isabl(kk)+iscalv(kk)+isablbord(kk)
202
203
204     enddo
205
206
207!!$
208!!$
209!!$            IF(H(i,j).gt.2.)  THEN
210!!$             if (mk(i,j).eq.1) then
211!!$               INF(2) = INF(2) + 1
212!!$               ISVOLF(2) = ISVOLF(2) + H(I,J)
213!!$             else
214!!$               INP(2) = INP(2) + 1
215!!$               ISVOL(2) = ISVOL(2) + H(I,J)
216!!$               ISACC(2) = ISACC(2) + ACC(I,J)
217!!$               ISABL(2) = ISABL(2) + ABL(I,J)
218!!$               ITJJA(2) = ITJJA(2) + TJULY(I,J)
219!!$               if (H(I,J).gt.HMAX_(2)) HMAX_(2)=H(I,J)
220!!$             endif
221!!$            ENDIF
222!!$            TACC(2) = TACC(2) + ACC(I,J)
223!!$             ISABLBORD(2) = ISABLBORD(2) + ABLBORD(I,J)
224!!$             ABLATOT(2) = ISABL(2) + ISABLBORD(2)
225!!$             ISCALV(2) = ISCALV(2) + CALV(I,J)
226!!$             ISBM(2) = ISACC(2)+ISABL(2)+ISCALV(2)+ISABLBORD(2)
227!!$             TBM(2) = TACC(2)+ISABL(2)+ISCALV(2)+ISABLBORD(2)
228!!$        ENDIF
229!!$
230!!$
231!!$      END DO 
232
233
234 
235do K=1,7
236
237! == Les moyennes     
238      IF(INP(K).ne.0)THEN
239          HMEAN_(K) = ISVOL(K)   /INP(K) ! /(DX*DY*INP(K))
240          BMEAN_(K) = ISBM(K)    /INP(K)
241          ACCMEAN_(K) = ISACC(K) /INP(K)
242          ABLMEAN_(K) = ISABL(K) /INP(K)
243          ABLBORDMEAN_(K) = ISABLBORD(K) /INP(K)
244          CALVMEAN_(K)    = ISCALV(K)    /INP(K)
245          ABLATOTMEAN_(K) = ABLATOT(K)   /INP(K)
246          ITJJAMEAN_(K)   = ITJJA(K)     /INP(K)
247       ENDIF
248! == Les volmes intergrées (3D)
249      ISVOL(K) = ISVOL(K)*DX*DY     
250              isacc(k)=isacc(k)*dx*dy
251              isabl(k)=isabl(k)*dx*dy
252              isablbord(k)=isablbord(k)*dx*dy
253              ablatot(k)=ablatot(k)*dx*dy
254              iscalv(k)=iscalv(k)*dx*dy
255              isbm(k)=isbm(k)*dx*dy
256              tacc(k)=tacc(k)*dx*dy
257              tbm(k)=tbm(k)*dx*dy
258enddo
259!nom_table='abl'
260!call printtable_r(abl,nom_table)
261
262!====================================================
263! -- Hemin40                : : 1   
264
265      if (NP.ne.0) then 
266!        HMEAN=VOL/NP
267!        VOL=VOL*DX*DY
268!        ISVOL(1) = VOL
269!        BMEAN=BMEAN/NP
270!        ACCMEAN=ACCMEAN/NP
271!        ITJJAMEAN = ITJJA(1)/NP
272!        ABLMEAN=BMEAN-ACCMEAN
273!!!     ABLMEAN = ABLMEAN/NP
274!        CALVMEAN=CALVMEAN/NP
275!!!     BMELTMEAN=BMELTMEAN/NP
276!!!     ABLBORDMEAN=ABLBORDMEAN/NP
277!!!     ABLAMEAN = ABLAMEAN/NP
278        TBMEAN=TBMEAN/NP
279        TBDOTMEAN=TBDOTMEAN/NP
280        VSMEAN=VSMEAN/NP
281!       VSDOTMEAN=VSDOTMEAN/NP
282        UZSMEAN=UZSMEAN/NP
283        UZSDOTMEAN=UZSDOTMEAN/NP
284        UZKMEAN=UZKMEAN/NP
285        DELTAVOL = HDOTMEAN*DX*DY
286        HDOTMEAN=HDOTMEAN/NP 
287    ELSE
288        TBMEAN=0.
289        TBDOTMEAN=0.
290        VSMEAN=0.
291        UZSMEAN=0.
292        UZSDOTMEAN=0.
293        UZKMEAN=0.
294        DELTAVOL = 0.
295        HDOTMEAN=0.
296      endif
297
298      BDOTMEAN=BDOTMEAN/NX/NY 
299
300! 2_writing outputs
301!------------------
302
303!         WRITE(35,904)TIME,VOL,DELTAVOL,NP,ISACC(1),ISABL(1),      &
304!          ISABLBORD(1),ABLATOT(1),ISCALV(1),ISBM(1),               &
305!          ISVOL(2),INPL,ISACC(2),                                  &
306!          ISABL(2),ISABLBORD(2),ABLATOT(2),ISCALV(2),ISBM(2),      &
307!          ISVOL(3),INPG,ISACC(3),ISABL(3),ISABLBORD(3),            &
308!          ABLATOT(3),ISCALV(3),ISBM(3),ISVOL(4),INPF,ISACC(4),     &
309!          ISABL(4),ISABLBORD(4),ABLATOT(4),ISCALV(4),ISBM(4),      &
310!          nint(HMEAN),nint(HMAX),BMEAN,TBMEAN,nint(VSMEAN),        &
311!          TBDOTMEAN,HDOTMEAN,BDOTMEAN,BMELTMEAN,NPAB,    &
312!!          TBDOTMEAN,VSDOTMEAN,HDOTMEAN,BDOTMEAN,BMELTMEAN,NPAB,    &
313!          NPCAL,HDOTMEAN_G,TACC(1),TBM(1),TACC(2),TBM(2),          &
314!          TACC(3),TBM(3),TACC(4),TBM(4),                           &
315!           ITJJAMEAN,ITJJAMEAN_L,ITJJAMEAN_G,ITJJAMEAN_F,          &
316!           ISVOL(5),INPLAB,ISACC(5),ABLATOT(5),ITJJAMEAN_LAB,      &
317!           ISVOL(6),INPKEEW,ISACC(6),ABLATOT(6),ITJJAMEAN_KEEW,    &
318!           ISVOL(7),INPINN,ISACC(7),ABLATOT(7),ITJJAMEAN_INN,      &
319!           ISVOL(8),INPCORD,ISACC(8),ABLATOT(8),ITJJAMEAN_CORD,    &
320!           ISVOL(9),INPSCAN,ISACC(9),ABLATOT(9),ITJJAMEAN_SCAN,    &
321!           ISVOL(10),INPBAR,ISACC(10),ABLATOT(10),ITJJAMEAN_BAR,   &
322!           ISVOL(11),INPKARA,ISACC(11),ABLATOT(11),ITJJAMEAN_KARA, &     
323!           ISVOL(12),INPSIB,ISACC(12),ABLATOT(12),ITJJAMEAN_SIB,   &     
324!           ISVOL(13),INPBRIT,ISACC(13),ABLATOT(13),ITJJAMEAN_BRIT, &
325!           HMEAN,HMEAN_L,HMEAN_G,HMEAN_F,HMEAN_LAB,HMEAN_KEEW,     &
326!           HMEAN_INN,HMEAN_CORD,HMEAN_SCAN,HMEAN_BAR,HMEAN_KARA,   &
327!           HMEAN_SIB,HMEAN_BRIT,                                   &
328!           HMAX,HMAX_L,HMAX_G,HMAX_F,HMAX_LAB,HMAX_KEEW,HMAX_INN,  &
329!           HMAX_CORD, HMAX_SCAN,HMAX_BAR,HMAX_KARA,HMAX_SIB,HMAX_BRIT     
330
331
332
333
334!print*,'========================='
335!print*,'ecrit short, time = ',time,dt
336!print*,'hdot ISABLBORD',hdot(30:40,100),ISABLBORD(2)
337
338          write(num_ritz,905)   time,isvol(1),deltavol,inp(1),isvolf(1),inf(1),isacc(1),isabl(1),    &   ! 8
339              ISABLBORD(1),ABLATOT(1),ISCALV(1),ISBM(1),     &                                     ! 4
340  isvol(2),inp(2),isvolf(2),inf(2),isacc(2),isabl(2),isablbord(2),ablatot(2),iscalv(2),isbm(2),itjjamean_(2),hmean_(2),Hmax_(2), & ! 13
341  isvol(3),inp(3),isvolf(3),inf(3),isacc(3),isabl(3),isablbord(3),ablatot(3),iscalv(3),isbm(3),itjjamean_(3),hmean_(3),hmax_(3), & ! 13
342  isvol(4),inp(4),isvolf(4),inf(4),isacc(4),isabl(4),isablbord(4),ablatot(4),iscalv(4),isbm(4),itjjamean_(4),hmean_(4),hmax_(4), & ! 13
343  isvol(5),inp(5),isvolf(5),inf(5),isacc(5),isabl(5),isablbord(5),ablatot(5),iscalv(5),isbm(5),itjjamean_(5),hmean_(5),hmax_(5), & ! 13
344  isvol(6),inp(6),isvolf(6),inf(6),isacc(6),isabl(6),isablbord(6),ablatot(6),iscalv(6),isbm(6),itjjamean_(6),hmean_(6),hmax_(6), & ! 13
345  isvol(7),inp(7),isvolf(7),inf(7),isacc(7),isabl(7),isablbord(7),ablatot(7),iscalv(7),isbm(7),itjjamean_(7),hmean_(7),hmax_(7), & ! 13
346  diff_H, water_bilan, sum(calv_dtt(:,:))/dtt, sum(ablbord_dtt(:,:))/dtt, sum(Bm_dtt(:,:))/dtt, sum(bmelt_dtt(:,:))/dtt !6
347
348 
349! !! format 900 faux, a reprendre
350
351905 format(f9.1,1x, e11.4,1x,e11.5, 1x,i5, 1x,e10.4,1x,i5 , 6(1x,e12.5), &
352         6(2(1x,e10.4,1x,i5), 9(1x,e11.4) ) , 6(1x,e15.8)) 
353 
354   
355940   format('%%%% ',a,'   time=',f8.0,' %%%%')
356end subroutine shortoutput
357
358
359end module output_laure16_mod
Note: See TracBrowser for help on using the repository browser.