GCC Code Coverage Report


Directory: ./
File: phys/simu_airs.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 485 0.0%
Branches: 0 666 0.0%

Line Branch Exec Source
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
1307