source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/simu_airs.F90 @ 5451

Last change on this file since 5451 was 2585, checked in by idelkadi, 9 years ago

Corrections des erreurs pour l'implementation du simulateur AIRS

File size: 35.3 KB
Line 
1       
2        module m_simu_airs
3
4        implicit none
5
6        real, parameter :: tau_thresh = 0.05 ! seuil nuages detectables
7        real, parameter :: p_thresh = 445.   ! seuil nuages hauts
8        real, parameter :: em_min = 0.2      ! seuils nuages semi-transp
9        real, parameter :: em_max = 0.85
10        real, parameter :: undef = 999.
11
12        contains
13
14        real function search_tropopause(P,T,alt,N) result(P_tropo)
15! this function searches for the tropopause pressure in [hPa].
16! The search is based on ideology described in
17! Reichler et al., Determining the tropopause height from gridded data,
18! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003
19
20        integer N,i,i_lev,first_point,exit_flag,i_dir
21        real P(N),T(N),alt(N),slope(N)
22        real P_min, P_max, slope_limit,slope_2km, &
23     & delta_alt_limit,tmp,delta_alt
24        parameter(P_min=75.0, P_max=470.0)   ! hPa
25        parameter(slope_limit=0.002)         ! 2 K/km converted to K/m
26        parameter(delta_alt_limit=2000.0)    ! 2000 meters
27
28        do i=1,N-1
29        slope(i)=-(T(i+1)-T(i))/(alt(i+1)-alt(i))
30        end do
31        slope(N)=slope(N-1)
32
33        if (P(1).gt.P(N)) then
34        i_dir= 1
35        i=1
36        else
37        i_dir=-1
38        i=N
39        end if
40
41        first_point=0
42        exit_flag=0
43        do while(exit_flag.eq.0)
44        if (P(i).ge.P_min.and.P(i).le.P_max) then
45        if (first_point.gt.0) then
46        delta_alt=alt(i)-alt(first_point)
47        if (delta_alt.ge.delta_alt_limit) then
48        slope_2km=(T(first_point)-T(i))/delta_alt
49        if (slope_2km.lt.slope_limit) then
50        exit_flag=1
51        else
52        first_point=0
53        end if
54        end if
55        end if
56        if (first_point.eq.0.and.slope(i).lt.slope_limit) first_point=i
57        end if
58        i=i+i_dir
59        if (i.le.1.or.i.ge.N) exit_flag=1
60        end do
61
62        if (first_point.le.0) P_tropo=65.4321
63        if (first_point.eq.1) P_tropo=64.5432
64        if (first_point.gt.1) then
65        tmp=(slope_limit-slope(first_point))/(slope(first_point+1)- &
66     & slope(first_point))*(P(first_point+1)-P(first_point))
67        P_tropo=P(first_point)+tmp
68        ! print*, 'P_tropo= ', tmp, P(first_point), P_tropo
69        end if
70
71! Ajout Marine
72        if (P_tropo .lt. 60. .or. P_tropo .gt. 470.) then
73        P_tropo = 999.
74        endif
75
76        return
77        end function search_tropopause
78
79
80
81        subroutine cloud_structure(len_cs, rneb_cs, temp_cs, &
82     & emis_cs, iwco_cs, &
83     & pres_cs, dz_cs, rhodz_cs, rad_cs, &
84     & cc_tot_cs, cc_hc_cs, cc_hist_cs, &
85     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
86     & pcld_hc_cs, tcld_hc_cs, &
87     & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
88     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
89     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
90     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
91     & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
92
93
94     
95        integer :: i, n, nss
96
97        integer, intent(in) :: len_cs
98        real, dimension(:), intent(in) :: rneb_cs, temp_cs
99        real, dimension(:), intent(in) :: emis_cs, iwco_cs, rad_cs
100        real, dimension(:), intent(in) :: pres_cs, dz_cs, rhodz_cs
101
102        real, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
103     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
104     & pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, &
105     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
106     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
107     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
108     & em_hist_cs, iwp_hist_cs, &
109     & deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs
110
111        real, dimension(len_cs) :: rneb_ord
112        real :: rneb_min
113
114        real, dimension(:), allocatable :: s, s_hc, s_hist, rneb_max
115        real, dimension(:), allocatable :: sCb, sThCi, sAnv
116        real, dimension(:), allocatable :: iwp_ss, pcld_ss, tcld_ss,&
117     & emis_ss
118        real, dimension(:), allocatable :: deltaz_ss, rad_ss
119
120
121        write(*,*) 'dans cloud_structure'
122
123        call ordonne(len_cs, rneb_cs, rneb_ord)
124       
125
126! Definition des sous_sections
127
128        rneb_min = rneb_ord(1)
129        nss = 1
130
131        do i = 2, size(rneb_cs)
132        if (rneb_ord(i) .gt. rneb_min) then
133        nss = nss+1
134        rneb_min = rneb_ord(i)
135        endif
136        enddo
137
138        allocate (s(nss))
139        allocate (s_hc(nss))
140        allocate (s_hist(nss))
141        allocate (rneb_max(nss))
142        allocate (emis_ss(nss))
143        allocate (pcld_ss(nss))
144        allocate (tcld_ss(nss))
145        allocate (iwp_ss(nss))
146        allocate (deltaz_ss(nss))
147        allocate (rad_ss(nss))
148        allocate (sCb(nss))
149        allocate (sThCi(nss))
150        allocate (sAnv(nss))
151
152        rneb_min = rneb_ord(1)
153        n = 1
154        s(1) = rneb_ord(1)
155        s_hc(1) = rneb_ord(1)
156        s_hist(1) = rneb_ord(1)
157        sCb(1) = rneb_ord(1)
158        sThCi(1) = rneb_ord(1)
159        sAnv(1) = rneb_ord(1)
160
161        rneb_max(1) = rneb_ord(1)
162
163        do i = 2, size(rneb_cs)
164        if (rneb_ord(i) .gt. rneb_min) then
165        n = n+1
166        s(n) = rneb_ord(i)-rneb_min
167        s_hc(n) = rneb_ord(i)-rneb_min
168        s_hist(n) = rneb_ord(i)-rneb_min
169        sCb(n) = rneb_ord(i)-rneb_min
170        sThCi(n) = rneb_ord(i)-rneb_min
171        sAnv(n) = rneb_ord(i)-rneb_min
172
173        rneb_max(n) = rneb_ord(i)
174        rneb_min = rneb_ord(i)
175        endif
176        enddo
177
178! Appel de sous_section
179
180        do i = 1, nss
181         call sous_section(len_cs, rneb_cs, temp_cs, &
182     &  emis_cs, iwco_cs, &
183     &  pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
184     &  rneb_max(i),s(i),s_hc(i),s_hist(i), &
185     &  sCb(i), sThCi(i), sAnv(i), &
186     &  emis_ss(i), &
187     &  pcld_ss(i), tcld_ss(i), iwp_ss(i), deltaz_ss(i), rad_ss(i))
188        enddo
189
190! Caracteristiques de la structure nuageuse
191
192        cc_tot_cs = 0.
193        cc_hc_cs = 0.
194        cc_hist_cs = 0.
195
196        cc_Cb_cs = 0.
197        cc_ThCi_cs = 0.
198        cc_Anv_cs = 0.
199
200        em_hc_cs = 0.
201        iwp_hc_cs = 0.
202        deltaz_hc_cs = 0.
203
204        em_hist_cs = 0.
205        iwp_hist_cs = 0.
206        deltaz_hist_cs = 0.
207        rad_hist_cs = 0.
208
209        pcld_hc_cs = 0.
210        tcld_hc_cs = 0.
211
212        pcld_Cb_cs = 0.
213        tcld_Cb_cs = 0.
214        em_Cb_cs = 0.
215
216        pcld_ThCi_cs = 0.
217        tcld_ThCi_cs = 0.
218        em_ThCi_cs = 0.
219
220        pcld_Anv_cs = 0.
221        tcld_Anv_cs = 0.
222        em_Anv_cs = 0.
223
224         do i = 1, nss
225
226        cc_tot_cs = cc_tot_cs + s(i)
227        cc_hc_cs = cc_hc_cs + s_hc(i)
228        cc_hist_cs = cc_hist_cs + s_hist(i)
229
230        cc_Cb_cs = cc_Cb_cs + sCb(i)
231        cc_ThCi_cs = cc_ThCi_cs + sThCi(i)
232        cc_Anv_cs = cc_Anv_cs + sAnv(i)
233
234        iwp_hc_cs = iwp_hc_cs + s_hc(i)*iwp_ss(i)
235        em_hc_cs = em_hc_cs + s_hc(i)*emis_ss(i)
236        deltaz_hc_cs = deltaz_hc_cs + s_hc(i)*deltaz_ss(i)
237
238        iwp_hist_cs = iwp_hist_cs + s_hist(i)*iwp_ss(i)
239        em_hist_cs = em_hist_cs + s_hist(i)*emis_ss(i)
240        deltaz_hist_cs = deltaz_hist_cs + s_hist(i)*deltaz_ss(i)
241        rad_hist_cs = rad_hist_cs + s_hist(i)*rad_ss(i)
242
243        pcld_hc_cs = pcld_hc_cs + s_hc(i)*pcld_ss(i)
244        tcld_hc_cs = tcld_hc_cs + s_hc(i)*tcld_ss(i)
245
246        pcld_Cb_cs = pcld_Cb_cs + sCb(i)*pcld_ss(i)
247        tcld_Cb_cs = tcld_Cb_cs + sCb(i)*tcld_ss(i)
248        em_Cb_cs = em_Cb_cs + sCb(i)*emis_ss(i)
249
250        pcld_ThCi_cs = pcld_ThCi_cs + sThCi(i)*pcld_ss(i)
251        tcld_ThCi_cs = tcld_ThCi_cs + sThCi(i)*tcld_ss(i)
252        em_ThCi_cs = em_ThCi_cs + sThCi(i)*emis_ss(i)
253
254        pcld_Anv_cs = pcld_Anv_cs + sAnv(i)*pcld_ss(i)
255        tcld_Anv_cs = tcld_Anv_cs + sAnv(i)*tcld_ss(i)
256        em_Anv_cs = em_Anv_cs + sAnv(i)*emis_ss(i)
257
258        enddo
259
260        deallocate(s)
261        deallocate (s_hc)
262        deallocate (s_hist)
263        deallocate (rneb_max)
264        deallocate (emis_ss)
265        deallocate (pcld_ss)
266        deallocate (tcld_ss)
267        deallocate (iwp_ss)
268        deallocate (deltaz_ss)
269        deallocate (rad_ss)
270        deallocate (sCb)
271        deallocate (sThCi)
272        deallocate (sAnv)
273
274       call normal_undef(pcld_hc_cs,cc_hc_cs)
275       call normal_undef(tcld_hc_cs,cc_hc_cs)
276       call normal_undef(iwp_hc_cs,cc_hc_cs)
277       call normal_undef(em_hc_cs,cc_hc_cs)
278       call normal_undef(deltaz_hc_cs,cc_hc_cs)
279       
280       call normal_undef(iwp_hist_cs,cc_hist_cs)
281       call normal_undef(em_hist_cs,cc_hist_cs)
282       call normal_undef(deltaz_hist_cs,cc_hist_cs)
283       call normal_undef(rad_hist_cs,cc_hist_cs)
284
285       call normal_undef(pcld_Cb_cs,cc_Cb_cs)
286       call normal_undef(tcld_Cb_cs,cc_Cb_cs)
287       call normal_undef(em_Cb_cs,cc_Cb_cs)
288
289       call normal_undef(pcld_ThCi_cs,cc_ThCi_cs)
290       call normal_undef(tcld_ThCi_cs,cc_ThCi_cs)
291       call normal_undef(em_ThCi_cs,cc_ThCi_cs)
292   
293       call normal_undef(pcld_Anv_cs,cc_Anv_cs)
294       call normal_undef(tcld_Anv_cs,cc_Anv_cs)
295       call normal_undef(em_Anv_cs,cc_Anv_cs)
296
297
298! Tests
299
300        if (cc_tot_cs .gt. maxval(rneb_cs) .and. &
301     & abs(cc_tot_cs-maxval(rneb_cs)) .gt. 1.e-4 )  then
302        write(*,*) 'cc_tot_cs > max rneb_cs'
303        write(*,*) cc_tot_cs, maxval(rneb_cs)
304        STOP
305        endif
306
307        if (iwp_hc_cs .lt. 0.) then
308        write(*,*) 'cloud_structure:'
309        write(*,*) 'iwp_hc_cs < 0'
310        STOP
311        endif
312 
313        end subroutine cloud_structure
314
315
316        subroutine normal_undef(num, den)
317
318        real, intent(in) :: den
319        real, intent(inout) :: num
320
321        if (den .ne. 0) then
322        num = num/den
323        else
324        num = undef
325        endif
326
327        end subroutine normal_undef
328
329
330        subroutine normal2_undef(res,num,den)
331
332        real, intent(in) :: den
333        real, intent(in) :: num
334        real, intent(out) :: res
335
336        if (den .ne. 0.) then
337        res = num/den
338        else
339        res = undef
340        endif
341
342        end subroutine normal2_undef
343
344
345        subroutine sous_section(len_cs, rneb_cs, temp_cs, &
346     & emis_cs, iwco_cs, &
347     & pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
348     & rnebmax, stot, shc, shist, &
349     & sCb, sThCi, sAnv, &
350     & emis, pcld, tcld, iwp, deltaz, rad)
351
352        integer, intent(in) :: len_cs
353        real, dimension(len_cs), intent(in) :: rneb_cs, temp_cs
354        real, dimension(len_cs), intent(in) :: emis_cs, iwco_cs, &
355     & rneb_ord
356        real, dimension(len_cs), intent(in) :: pres_cs, dz_cs, rad_cs
357        real, dimension(len_cs), intent(in) :: rhodz_cs
358        real, dimension(len_cs) :: tau_cs, w
359        real, intent(in) :: rnebmax
360        real, intent(inout) :: stot, shc, shist
361        real, intent(inout) :: sCb, sThCi, sAnv
362        real, intent(out) :: emis, pcld, tcld, iwp, deltaz, rad
363
364        integer :: i, ideb, ibeg, iend, nuage, visible
365        real :: som, som_tau, som_iwc, som_dz, som_rad, tau
366
367
368! Ponderation: 1 pour les nuages, 0 pour les trous
369
370        do i = 1, len_cs
371        if (rneb_cs(i) .ge. rnebmax) then
372        w(i) = 1.
373        else
374        w(i) = 0.
375        endif
376        enddo
377
378! Calcul des epaisseurs optiques a partir des emissivites
379
380        som = 0.
381        do i = 1, len_cs
382        if (emis_cs(i) .eq. 1.) then
383        tau_cs(i) = 10.
384        else
385        tau_cs(i) = -log(1.-emis_cs(i))
386        endif
387        som = som+tau_cs(i)
388        enddo
389
390
391        ideb = 1
392        nuage = 0
393        visible = 0
394
395
396! Boucle sur les nuages
397        do while (ideb .ne. 0 .and. ideb .le. len_cs)   
398
399
400! Definition d'un nuage de la sous-section
401
402        call topbot(ideb, w, ibeg, iend)
403        ideb = iend+1
404
405        if (ibeg .gt. 0) then
406
407        nuage = nuage + 1
408
409! On determine les caracteristiques du nuage
410! (ep. optique, ice water path, pression, temperature)
411
412        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
413     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
414     & som_tau, som_iwc, som_dz, som_rad)
415
416! On masque le nuage s'il n'est pas detectable
417
418        call masque(ibeg, iend, som_tau, visible, w)
419
420        endif
421
422! Fin boucle nuages
423        enddo
424
425
426! Analyse du nuage detecte
427
428        call topbot(1, w, ibeg, iend)
429
430        if (ibeg .gt. 0) then
431
432        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
433     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
434     & som_tau, som_iwc, som_dz, som_rad)
435
436        tau = som_tau
437        emis = 1. - exp(-tau)
438        iwp = som_iwc
439        deltaz = som_dz
440        rad = som_rad
441
442        if (pcld .gt. p_thresh) then
443
444        shc = 0.
445        shist = 0.
446        sCb = 0.
447        sThCi = 0.
448        sAnv = 0.
449
450        else
451
452        if (emis .lt. em_min .or. emis .gt. em_max  &
453     & .or. tcld .gt. 230.) then
454        shist = 0.
455        endif
456
457        if (emis .lt. 0.98) then
458        sCb = 0.
459        endif
460
461        if (emis .gt. 0.5 .or. emis .lt. 0.1) then
462        sThCi = 0.
463        endif
464
465        if (emis .le. 0.5 .or. emis .ge. 0.98) then
466        sAnv = 0.
467        endif
468
469        endif
470
471        else
472
473        tau = 0.
474        emis = 0.
475        iwp = 0.
476        deltaz = 0.
477        pcld = 0.
478        tcld = 0.
479        stot = 0.
480        shc = 0.
481        shist = 0.
482        rad = 0.
483        sCb = 0.
484        sThCi = 0.
485        sAnv = 0.
486
487        endif
488
489
490! Tests
491
492        if (iwp .lt. 0.) then
493        write(*,*) 'ideb iwp =', ideb, iwp
494        STOP
495        endif
496
497        if (deltaz .lt. 0.) then
498        write(*,*) 'ideb deltaz =', ideb, deltaz
499        STOP
500        endif
501
502        if (emis .lt. 0.048 .and. emis .ne. 0.) then
503        write(*,*) 'ideb emis =', ideb, emis
504        STOP
505        endif
506
507        end subroutine sous_section
508
509
510        subroutine masque (ibeg, iend, som_tau, &
511     & visible, w)
512
513        integer, intent(in) :: ibeg, iend
514        real, intent(in) :: som_tau
515
516        integer, intent(inout) :: visible
517        real, dimension(:), intent(inout) :: w
518
519        integer :: i
520
521
522
523! Masque
524
525! Cas ou il n'y a pas de nuage visible au-dessus
526
527        if (visible .eq. 0) then
528
529        if (som_tau .lt. tau_thresh) then
530        do i = ibeg, iend
531        w(i) = 0.
532        enddo
533        else
534        visible = 1
535        endif
536
537! Cas ou il y a un nuage visible au-dessus
538
539        else
540
541        do i = ibeg, iend
542        w(i) = 0.
543        enddo
544
545        endif
546
547
548        end subroutine masque
549
550
551         subroutine caract (ibeg, iend, temp_cs, tau_cs, iwco_cs, &
552     & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
553     & som_tau, som_iwc, som_dz, som_rad)
554
555        integer, intent(in) :: ibeg, iend
556        real, dimension(:), intent(in) :: tau_cs, iwco_cs, temp_cs
557        real, dimension(:), intent(in) :: pres_cs, dz_cs, rad_cs
558        real, dimension(:), intent(in) :: rhodz_cs
559        real, intent(out) :: som_tau, som_iwc, som_dz, som_rad
560        real , intent(out) :: pcld, tcld
561
562        integer :: i, ibase, imid
563
564! Somme des epaisseurs optiques et des contenus en glace sur le nuage
565
566        som_tau = 0.
567        som_iwc = 0.
568        som_dz = 0.
569        som_rad = 0.
570        ibase = -100
571
572        do i = ibeg, iend
573
574        som_tau = som_tau + tau_cs(i)
575
576        som_dz = som_dz + dz_cs(i)
577        som_iwc = som_iwc + iwco_cs(i)*1000*rhodz_cs(i)  ! en g/m2
578        som_rad = som_rad + rad_cs(i)*dz_cs(i)
579
580        if (som_tau .gt. 3. .and. ibase .eq. -100) then
581        ibase = i-1
582        endif
583
584        enddo
585
586        if (som_dz .ne. 0.) then
587        som_rad = som_rad/som_dz
588        else
589        write(*,*) 'som_dez = 0 STOP'
590        write(*,*) 'ibeg, iend =', ibeg, iend
591        do i = ibeg, iend
592        write(*,*) dz_cs(i), rhodz_cs(i)
593        enddo
594        STOP
595        endif
596
597! Determination de Pcld et Tcld
598
599       if (ibase .lt. ibeg) then
600       ibase = ibeg
601       endif
602
603       imid = (ibeg+ibase)/2
604
605       pcld = pres_cs(imid)/100.        ! pcld en hPa
606       tcld = temp_cs(imid)
607
608
609       end subroutine caract
610 
611        subroutine topbot(ideb,w,ibeg,iend)
612
613        integer, intent(in) :: ideb
614        real, dimension(:), intent(in) :: w
615        integer, intent(out) :: ibeg, iend
616
617        integer :: i, itest
618
619        itest = 0
620        ibeg = 0
621        iend = 0
622
623        do i = ideb, size(w)
624
625        if (w(i) .eq. 1. .and. itest .eq. 0) then
626        ibeg = i
627        itest = 1
628        endif
629
630        enddo
631
632
633        i = ibeg
634        do while (w(i) .eq. 1. .and. i .le. size(w))
635        i = i+1
636        enddo
637        iend = i-1
638
639
640        end subroutine topbot
641
642        subroutine ordonne(len_cs, rneb_cs, rneb_ord)
643
644        integer, intent(in) :: len_cs
645        real, dimension(:), intent(in) :: rneb_cs
646        real, dimension(:), intent(out) :: rneb_ord
647
648        integer :: i, j, ind_min
649
650        real, dimension(len_cs) :: rneb
651        real :: rneb_max
652
653
654        do i = 1, size(rneb_cs)
655        rneb(i) = rneb_cs(i)
656        enddo
657
658        do j = 1, size(rneb_cs)
659
660        rneb_max = 100.
661
662        do i = 1, size(rneb_cs)
663        if (rneb(i) .lt. rneb_max) then
664        rneb_max = rneb(i)
665        ind_min = i
666        endif
667        enddo
668
669        rneb_ord(j) = rneb_max
670        rneb(ind_min) = 100.
671
672        enddo
673       
674        end subroutine ordonne
675
676 
677        subroutine sim_mesh(rneb_1D, temp_1D, emis_1D, &
678     & iwcon_1D, rad_1D, &
679     & pres, dz, &
680     & rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh,&
681     & tcld_hc_mesh, &
682     & em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
683     & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
684     & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
685     & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
686     & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
687     & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
688
689       USE dimphy
690
691       real, dimension(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
692     & iwcon_1D, rad_1D
693        real, dimension(klev), intent(in) :: pres, dz, rhodz_1D
694        real, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
695        real, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
696        real, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
697     & iwp_hc_mesh
698
699        real, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
700        real, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
701     & em_ThCi_mesh
702        real, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
703
704        real, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
705        real, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh
706
707        real, dimension(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
708     & iwco_cs
709        real, dimension(:), allocatable :: pres_cs, dz_cs, rad_cs, &
710     & rhodz_cs
711
712        integer :: i,j,l
713        integer :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics
714
715        real :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,&
716     & som_hist
717        real :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
718     & som_deltaz_hist
719        real :: som_rad_hist
720        real :: som_Cb, som_ThCi, som_Anv
721        real :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
722        real :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
723        real :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
724        real :: tsom_tot, tsom_hc, tsom_hist
725        real :: prod, prod_hh
726
727        real :: cc_tot_cs, cc_hc_cs, cc_hist_cs
728        real :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
729        real :: pcld_hc_cs, tcld_hc_cs
730        real :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
731        real :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
732        real :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
733        real :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
734        real :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs
735
736        real, dimension(klev) :: test_tot, test_hc, test_hist
737        real, dimension(klev) :: test_pcld, test_tcld, test_em, test_iwp
738
739
740        do j = 1, klev
741        write(*,*) 'simu_airs, j, rneb_1D =', rneb_1D(j)
742        enddo
743
744! Definition des structures nuageuses, de la plus haute a la plus basse
745
746        num_cs = 0
747        ltop = klev-1
748
749        prod = 1.
750
751        som_emi_hc = 0.
752        som_emi_hist = 0.
753        som_pcl_hc = 0.
754        som_tcl_hc = 0.
755        som_iwp_hc = 0.
756        som_iwp_hist = 0.
757        som_deltaz_hc = 0.
758        som_deltaz_hist = 0.
759        som_rad_hist = 0.
760        som_hc = 0.
761        som_hist = 0.
762
763        som_Cb = 0.
764        som_ThCi = 0.
765        som_Anv = 0.
766
767        som_emi_Cb = 0.
768        som_pcld_Cb = 0.
769        som_tcld_Cb = 0.
770
771        som_emi_ThCi = 0.
772        som_pcld_ThCi = 0.
773        som_tcld_ThCi = 0.
774
775        som_emi_Anv = 0.
776        som_pcld_Anv = 0.
777        som_tcld_Anv = 0.
778
779        tsom_tot = 0.
780        tsom_hc = 0.
781        tsom_hist = 0.
782
783        do while (ltop .ge. 1)   ! Boucle sur toute la colonne
784
785        itop = 0
786
787        do l = ltop,1,-1
788
789        if (itop .eq. 0 .and. rneb_1D(l) .gt. 0.001 ) then
790        itop = l
791        endif
792
793        enddo
794
795        ibot = itop
796
797        do while (rneb_1D(ibot) .gt. 0.001 .and. ibot .ge. 1)
798        ibot = ibot -1
799        enddo
800
801
802        ibot = ibot+1
803
804        if (itop .gt. 0) then    ! itop > 0
805
806        num_cs = num_cs +1
807        len_cs = itop-ibot+1
808
809! Allocation et definition des variables de la structure nuageuse
810! le premier indice denote ce qui est le plus haut
811
812        allocate (rneb_cs(len_cs))
813        allocate (temp_cs(len_cs))
814        allocate (emis_cs(len_cs))
815        allocate (iwco_cs(len_cs))
816        allocate (pres_cs(len_cs))
817        allocate (dz_cs(len_cs))
818        allocate (rad_cs(len_cs))
819        allocate (rhodz_cs(len_cs))
820
821        ics = 0
822
823        do i = itop, ibot, -1
824        ics = ics + 1
825        rneb_cs(ics) = rneb_1D(i)
826        temp_cs(ics) = temp_1D(i)
827        emis_cs(ics) = emis_1D(i)
828        iwco_cs(ics) = iwcon_1D(i)
829        rad_cs(ics) = rad_1D(i)
830        pres_cs(ics) = pres(i)
831        dz_cs(ics) = dz(i)
832        rhodz_cs(ics) = rhodz_1D(i)
833        enddo
834
835! Appel du sous_programme cloud_structure
836
837        call cloud_structure(len_cs,rneb_cs,temp_cs,emis_cs,iwco_cs,&
838     & pres_cs, dz_cs, rhodz_cs, rad_cs, &
839     & cc_tot_cs, cc_hc_cs, cc_hist_cs, &
840     & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
841     & pcld_hc_cs, tcld_hc_cs, &
842     & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
843     & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
844     & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
845     & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
846     & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)
847
848
849        deallocate (rneb_cs)
850        deallocate (temp_cs)
851        deallocate (emis_cs)
852        deallocate (iwco_cs)
853        deallocate (pres_cs)
854        deallocate (dz_cs)
855        deallocate (rad_cs)
856        deallocate (rhodz_cs)
857
858
859! Pour la couverture nuageuse sur la maille
860
861        prod_hh = prod
862
863        prod = prod*(1.-cc_tot_cs)
864
865! Pour les autres variables definies sur la maille
866
867        som_emi_hc = som_emi_hc + em_hc_cs*cc_hc_cs*prod_hh
868        som_iwp_hc = som_iwp_hc + iwp_hc_cs*cc_hc_cs*prod_hh
869        som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs*cc_hc_cs*prod_hh
870
871        som_emi_Cb = som_emi_Cb + em_Cb_cs*cc_Cb_cs*prod_hh
872        som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs*cc_Cb_cs*prod_hh
873        som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs*cc_Cb_cs*prod_hh
874
875        som_emi_ThCi = som_emi_ThCi + em_ThCi_cs*cc_ThCi_cs*prod_hh
876        som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs*cc_ThCi_cs*prod_hh
877        som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs*cc_ThCi_cs*prod_hh
878
879        som_emi_Anv = som_emi_Anv + em_Anv_cs*cc_Anv_cs*prod_hh
880        som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs*cc_Anv_cs*prod_hh
881        som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs*cc_Anv_cs*prod_hh
882
883        som_emi_hist = som_emi_hist + em_hist_cs*cc_hist_cs*prod_hh
884        som_iwp_hist = som_iwp_hist + iwp_hist_cs*cc_hist_cs*prod_hh
885        som_deltaz_hist = som_deltaz_hist + &
886     & deltaz_hist_cs*cc_hist_cs*prod_hh
887        som_rad_hist = som_rad_hist + rad_hist_cs*cc_hist_cs*prod_hh
888
889        som_pcl_hc = som_pcl_hc + pcld_hc_cs*cc_hc_cs*prod_hh
890        som_tcl_hc = som_tcl_hc + tcld_hc_cs*cc_hc_cs*prod_hh
891
892        som_hc = som_hc + cc_hc_cs*prod_hh
893        som_hist = som_hist + cc_hist_cs*prod_hh
894
895        som_Cb = som_Cb + cc_Cb_cs*prod_hh
896        som_ThCi = som_ThCi + cc_ThCi_cs*prod_hh
897        som_Anv = som_Anv + cc_Anv_cs*prod_hh
898
899
900! Pour test
901       
902        call test_bornes('cc_tot_cs     ',cc_tot_cs,1.,0.)
903        call test_bornes('cc_hc_cs      ',cc_hc_cs,1.,0.)
904        call test_bornes('cc_hist_cs    ',cc_hist_cs,1.,0.)
905        call test_bornes('pcld_hc_cs    ',pcld_hc_cs,1200.,0.)
906        call test_bornes('tcld_hc_cs    ',tcld_hc_cs,1000.,100.)
907        call test_bornes('em_hc_cs      ',em_hc_cs,1000.,0.048)
908
909        test_tot(num_cs) = cc_tot_cs
910        test_hc(num_cs) = cc_hc_cs
911        test_hist(num_cs) = cc_hist_cs
912        test_pcld(num_cs) = pcld_hc_cs
913        test_tcld(num_cs) = tcld_hc_cs
914        test_em(num_cs) = em_hc_cs
915        test_iwp(num_cs) = iwp_hc_cs
916
917        tsom_tot = tsom_tot + cc_tot_cs
918        tsom_hc = tsom_hc + cc_hc_cs
919        tsom_hist = tsom_hist + cc_hist_cs
920
921
922        endif                   ! itop > 0
923
924        ltop = ibot -1
925
926        enddo                   ! fin de la boucle sur la colonne
927
928        N_CS = num_cs
929
930
931! Determination des variables de sortie
932
933        if (N_CS .gt. 0) then   ! if N_CS>0
934
935        cc_tot_mesh = 1. - prod
936
937        cc_hc_mesh = som_hc
938        cc_hist_mesh = som_hist
939
940        cc_Cb_mesh = som_Cb
941        cc_ThCi_mesh = som_ThCi
942        cc_Anv_mesh = som_Anv
943
944        call normal2_undef(pcld_hc_mesh,som_pcl_hc, &
945     & cc_hc_mesh)
946        call normal2_undef(tcld_hc_mesh,som_tcl_hc, &
947     & cc_hc_mesh)
948        call normal2_undef(em_hc_mesh,som_emi_hc, &
949     & cc_hc_mesh)
950        call normal2_undef(iwp_hc_mesh,som_iwp_hc, &
951     & cc_hc_mesh)
952        call normal2_undef(deltaz_hc_mesh,som_deltaz_hc, &
953     & cc_hc_mesh)
954
955        call normal2_undef(em_Cb_mesh,som_emi_Cb, &
956     & cc_Cb_mesh)
957        call normal2_undef(tcld_Cb_mesh,som_tcld_Cb, &
958     & cc_Cb_mesh)
959        call normal2_undef(pcld_Cb_mesh,som_pcld_Cb, &
960     & cc_Cb_mesh)
961
962        call normal2_undef(em_ThCi_mesh,som_emi_ThCi, &
963     & cc_ThCi_mesh)
964        call normal2_undef(tcld_ThCi_mesh,som_tcld_ThCi, &
965     & cc_ThCi_mesh)
966        call normal2_undef(pcld_ThCi_mesh,som_pcld_ThCi, &
967     & cc_ThCi_mesh)
968
969       call normal2_undef(em_Anv_mesh,som_emi_Anv, &
970     & cc_Anv_mesh)
971        call normal2_undef(tcld_Anv_mesh,som_tcld_Anv, &
972     & cc_Anv_mesh)
973        call normal2_undef(pcld_Anv_mesh,som_pcld_Anv, &
974     & cc_Anv_mesh)
975
976
977        call normal2_undef(em_hist_mesh,som_emi_hist, &
978     & cc_hist_mesh)
979        call normal2_undef(iwp_hist_mesh,som_iwp_hist, &
980     & cc_hist_mesh)
981        call normal2_undef(deltaz_hist_mesh,som_deltaz_hist, &
982     & cc_hist_mesh)
983        call normal2_undef(rad_hist_mesh,som_rad_hist, &
984     & cc_hist_mesh)
985
986
987! Tests
988
989        ! Tests
990
991       if (cc_tot_mesh .gt. tsom_tot .and. &
992     & abs(cc_tot_mesh-tsom_tot) .gt. 1.e-4) then
993        write(*,*) 'cc_tot_mesh > tsom_tot'
994        write(*,*) cc_tot_mesh, tsom_tot
995        STOP
996        endif
997
998        if (cc_tot_mesh .lt. maxval(test_tot(1:N_CS)) .and. &
999     & abs(cc_tot_mesh-maxval(test_tot(1:N_CS))) .gt. 1.e-4) then
1000        write(*,*) 'cc_tot_mesh < max'
1001        write(*,*) cc_tot_mesh, maxval(test_tot(1:N_CS))
1002        STOP
1003        endif
1004
1005        if (cc_hc_mesh .gt. tsom_hc .and. &
1006     & abs(cc_hc_mesh-tsom_hc) .gt. 1.e-4) then
1007        write(*,*) 'cc_hc_mesh > tsom_hc'
1008        write(*,*) cc_hc_mesh, tsom_hc
1009        STOP
1010        endif
1011
1012        if (cc_hc_mesh .lt. maxval(test_hc(1:N_CS)) .and. &
1013     & abs(cc_hc_mesh-maxval(test_hc(1:N_CS))) .gt. 1.e-4) then
1014        write(*,*) 'cc_hc_mesh < max'
1015        write(*,*) cc_hc_mesh, maxval(test_hc(1:N_CS))
1016        STOP
1017        endif
1018
1019        if (cc_hist_mesh .gt. tsom_hist .and. &
1020     & abs(cc_hist_mesh-tsom_hist) .gt. 1.e-4) then
1021        write(*,*) 'cc_hist_mesh > tsom_hist'
1022        write(*,*) cc_hist_mesh, tsom_hist
1023        STOP
1024        endif
1025
1026        if (cc_hist_mesh .lt. 0.) then
1027        write(*,*) 'cc_hist_mesh < 0'
1028        write(*,*) cc_hist_mesh
1029        STOP
1030        endif
1031
1032        if ((pcld_hc_mesh .gt. maxval(test_pcld(1:N_CS)) .or. &
1033     & pcld_hc_mesh .lt. minval(test_pcld(1:N_CS))) .and. &
1034     & abs(pcld_hc_mesh-maxval(test_pcld(1:N_CS))) .gt. 1. .and. &
1035     & maxval(test_pcld(1:N_CS)) .ne. 999. &
1036     & .and. minval(test_pcld(1:N_CS)) .ne. 999.) then
1037        write(*,*) 'pcld_hc_mesh est faux'
1038        write(*,*) pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
1039     & minval(test_pcld(1:N_CS))
1040        STOP
1041        endif
1042
1043       if ((tcld_hc_mesh .gt. maxval(test_tcld(1:N_CS)) .or. &
1044     & tcld_hc_mesh .lt. minval(test_tcld(1:N_CS))) .and. &
1045     & abs(tcld_hc_mesh-maxval(test_tcld(1:N_CS))) .gt. 0.1 .and. &
1046     & maxval(test_tcld(1:N_CS)) .ne. 999. &
1047     & .and. minval(test_tcld(1:N_CS)) .ne. 999.) then
1048        write(*,*) 'tcld_hc_mesh est faux'
1049        write(*,*) tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
1050     & minval(test_tcld(1:N_CS))
1051        endif
1052
1053        if ((em_hc_mesh .gt. maxval(test_em(1:N_CS)) .or. &
1054     & em_hc_mesh .lt. minval(test_em(1:N_CS))) .and. &
1055     & abs(em_hc_mesh-maxval(test_em(1:N_CS))) .gt. 1.e-4 .and. &
1056     & minval(test_em(1:N_CS)) .ne. 999. .and. &
1057     & maxval(test_em(1:N_CS)) .ne. 999. ) then
1058        write(*,*) 'em_hc_mesh est faux'
1059        write(*,*) em_hc_mesh, maxval(test_em(1:N_CS)), &
1060     & minval(test_em(1:N_CS))
1061        STOP
1062        endif
1063
1064        else               ! if N_CS>0
1065
1066        cc_tot_mesh = 0.
1067        cc_hc_mesh = 0.
1068        cc_hist_mesh = 0.
1069
1070        cc_Cb_mesh = 0.
1071        cc_ThCi_mesh = 0.
1072        cc_Anv_mesh = 0.
1073
1074        iwp_hc_mesh = undef
1075        deltaz_hc_mesh = undef
1076        em_hc_mesh = undef
1077        iwp_hist_mesh = undef
1078        deltaz_hist_mesh = undef
1079        rad_hist_mesh = undef
1080        em_hist_mesh = undef
1081        pcld_hc_mesh = undef
1082        tcld_hc_mesh = undef
1083
1084        pcld_Cb_mesh = undef
1085        tcld_Cb_mesh = undef
1086        em_Cb_mesh = undef
1087
1088        pcld_ThCi_mesh = undef
1089        tcld_ThCi_mesh = undef
1090        em_ThCi_mesh = undef
1091
1092        pcld_Anv_mesh = undef
1093        tcld_Anv_mesh = undef
1094        em_Anv_mesh = undef
1095
1096
1097        endif                  ! if N_CS>0
1098
1099        end subroutine sim_mesh
1100
1101        subroutine test_bornes(sx,x,bsup,binf)
1102
1103        real, intent(in) :: x, bsup, binf
1104        character*14, intent(in) :: sx
1105
1106        if (x .gt. bsup .or. x .lt. binf) then
1107        write(*,*) sx, 'est faux'
1108        write(*,*) sx, x
1109        STOP
1110        endif
1111 
1112        end subroutine test_bornes
1113
1114        end module m_simu_airs
1115
1116
1117        subroutine simu_airs &
1118     & (itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
1119     & geop_airs, pplay_airs, paprs_airs, &
1120     & map_prop_hc,map_prop_hist,&
1121     & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
1122     & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
1123     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
1124     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
1125     & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
1126     & map_ntot,map_hc,map_hist,&
1127     & map_Cb,map_ThCi,map_Anv,alt_tropo )
1128
1129        USE dimphy
1130        USE m_simu_airs
1131
1132        IMPLICIT NONE
1133
1134        include "YOMCST.h"
1135
1136        integer,intent(in) :: itap
1137
1138        real, dimension(klon,klev), intent(in) :: &
1139     & rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
1140     & rad_airs, geop_airs, pplay_airs, paprs_airs
1141
1142       real, dimension(klon,klev) :: &
1143     & rhodz_airs, rho_airs, iwcon_airs
1144
1145        real, dimension(klon),intent(out) :: alt_tropo
1146
1147        real, dimension(klev) :: rneb_1D, temp_1D, &
1148     & emis_1D, rad_1D, pres_1D, alt_1D, &
1149     & rhodz_1D, dz_1D, iwcon_1D
1150
1151        integer :: i, j
1152
1153        real :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1154        real :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
1155        real :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1156        real :: em_hist_mesh, iwp_hist_mesh
1157        real :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
1158        real :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
1159        real :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
1160        real :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh
1161
1162        real, dimension(klon),intent(out) :: map_prop_hc, map_prop_hist
1163        real, dimension(klon),intent(out) :: map_emis_hc, map_iwp_hc
1164        real, dimension(klon),intent(out) :: map_deltaz_hc, map_pcld_hc
1165        real, dimension(klon),intent(out) :: map_tcld_hc
1166        real, dimension(klon),intent(out) :: map_emis_Cb,map_pcld_Cb,map_tcld_Cb
1167        real, dimension(klon),intent(out) :: &
1168     & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi
1169        real, dimension(klon),intent(out) :: &
1170     & map_emis_Anv,map_pcld_Anv,map_tcld_Anv
1171        real, dimension(klon),intent(out) :: &
1172     & map_emis_hist,map_iwp_hist,map_deltaz_hist,&
1173     & map_rad_hist
1174       real, dimension(klon),intent(out) :: map_ntot,map_hc,map_hist
1175       real, dimension(klon),intent(out) :: map_Cb,map_ThCi,map_Anv
1176 
1177 
1178        write(*,*) 'simu_airs'
1179        write(*,*) 'itap, klon, klev', itap, klon, klev
1180        write(*,*) 'RG, RD =', RG, RD
1181
1182
1183! Definition des variables 1D
1184
1185        do i = 1, klon
1186        do j = 1, klev-1
1187        rhodz_airs(i,j) = &
1188     & (paprs_airs(i,j)-paprs_airs(i,j+1))/RG
1189        enddo
1190        rhodz_airs(i,klev) = 0.
1191        enddo
1192
1193        do i = 1, klon
1194        do j = 1,klev
1195        rho_airs(i,j) = &
1196     & pplay_airs(i,j)/(temp_airs(i,j)*RD)
1197
1198        if (rneb_airs(i,j) .gt. 0.001) then
1199        iwcon_airs(i,j) = iwcon0_airs(i,j)/rneb_airs(i,j)
1200        else
1201        iwcon_airs(i,j) = 0.
1202        endif
1203 
1204        enddo
1205        enddo
1206
1207!=============================================================================
1208
1209        do i = 1, klon  ! boucle sur les points de grille
1210
1211!=============================================================================
1212       
1213        do j = 1,klev
1214
1215        rneb_1D(j) = rneb_airs(i,j)
1216        temp_1D(j) = temp_airs(i,j)
1217        emis_1D(j) = cldemi_airs(i,j)
1218        iwcon_1D(j) = iwcon_airs(i,j)
1219        rad_1D(j) = rad_airs(i,j)
1220        pres_1D(j) = pplay_airs(i,j)
1221        alt_1D(j) = geop_airs(i,j)/RG
1222        rhodz_1D(j) = rhodz_airs(i,j)
1223        dz_1D(j) = rhodz_airs(i,j)/rho_airs(i,j)
1224
1225        enddo
1226
1227        alt_tropo(i) = &
1228     & search_tropopause(pres_1D/100.,temp_1D,alt_1D,klev)
1229
1230
1231! Appel du ss-programme sim_mesh
1232
1233!        if (itap .eq. 1 ) then
1234
1235        call sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
1236     & pres_1D, dz_1D, rhodz_1D, &
1237     & cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
1238     & pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
1239     & deltaz_hc_mesh,&
1240     & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
1241     & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
1242     & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
1243     & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
1244     & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)
1245
1246         write(*,*) '===================================='
1247         write(*,*) 'itap, i:', itap, i
1248         write(*,*) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
1249     & iwp_hc, em_hist, iwp_hist ='
1250         write(*,*) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
1251         write(*,*) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
1252         write(*,*)  em_hist_mesh, iwp_hist_mesh
1253
1254!        endif
1255
1256! Definition des variables a ecrire dans le fichier de sortie
1257
1258        call normal2_undef(map_prop_hc(i),cc_hc_mesh, &
1259     & cc_tot_mesh)
1260        call normal2_undef(map_prop_hist(i),cc_hist_mesh, &
1261     & cc_tot_mesh)
1262
1263       map_emis_hc(i) = em_hc_mesh
1264       map_iwp_hc(i) = iwp_hc_mesh
1265       map_deltaz_hc(i) = deltaz_hc_mesh
1266       map_pcld_hc(i) = pcld_hc_mesh
1267       map_tcld_hc(i) = tcld_hc_mesh
1268
1269       map_emis_Cb(i) = em_Cb_mesh
1270       map_pcld_Cb(i) = pcld_Cb_mesh
1271       map_tcld_Cb(i) = tcld_Cb_mesh
1272
1273       map_emis_ThCi(i) = em_ThCi_mesh
1274       map_pcld_ThCi(i) = pcld_ThCi_mesh
1275       map_tcld_ThCi(i) = tcld_ThCi_mesh
1276
1277       map_emis_Anv(i) = em_Anv_mesh
1278       map_pcld_Anv(i) = pcld_Anv_mesh
1279       map_tcld_Anv(i) = tcld_Anv_mesh
1280
1281       map_emis_hist(i) = em_hist_mesh
1282       map_iwp_hist(i) = iwp_hist_mesh
1283       map_deltaz_hist(i) = deltaz_hist_mesh
1284       map_rad_hist(i) = rad_hist_mesh
1285
1286       map_ntot(i) = cc_tot_mesh
1287       map_hc(i) = cc_hc_mesh
1288       map_hist(i) = cc_hist_mesh
1289
1290       map_Cb(i) = cc_Cb_mesh
1291       map_ThCi(i) = cc_ThCi_mesh
1292       map_Anv(i) = cc_Anv_mesh
1293
1294
1295        enddo         ! fin boucle sur les points de grille
1296
1297       
1298
1299        end subroutine simu_airs
1300
Note: See TracBrowser for help on using the repository browser.