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