source: LMDZ6/branches/LMDZ-COSP/libf/phylmd/simu_airs.F90 @ 3847

Last change on this file since 3847 was 3531, checked in by Laurent Fairhead, 6 years ago

Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86
Still some work to be done in phylmd subdirectories

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