module m_simu_airs USE print_control_mod, ONLY: prt_level,lunout implicit none REAL, PARAMETER :: tau_thresh = 0.05 ! seuil nuages detectables REAL, PARAMETER :: p_thresh = 445. ! seuil nuages hauts REAL, PARAMETER :: em_min = 0.2 ! seuils nuages semi-transp REAL, PARAMETER :: em_max = 0.85 REAL, PARAMETER :: undef = 999. contains REAL function search_tropopause(P,T,alt,N) result(P_tropo) ! this function searches for the tropopause pressure in [hPa]. ! The search is based on ideology described in ! Reichler et al., Determining the tropopause height from gridded data, ! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003 INTEGER N,i,i_lev,first_point,exit_flag,i_dir REAL P(N),T(N),alt(N),slope(N) REAL P_min, P_max, slope_limit,slope_2km, & & delta_alt_limit,tmp,delta_alt PARAMETER(P_min=75.0, P_max=470.0) ! hPa PARAMETER(slope_limit=0.002) ! 2 K/km converted to K/m PARAMETER(delta_alt_limit=2000.0) ! 2000 meters do i=1,N-1 slope(i)=-(T(i+1)-T(i))/(alt(i+1)-alt(i)) end do slope(N)=slope(N-1) if (P(1).gt.P(N)) then i_dir= 1 i=1 else i_dir=-1 i=N end if first_point=0 exit_flag=0 do while(exit_flag.eq.0) if (P(i).ge.P_min.and.P(i).le.P_max) then if (first_point.gt.0) then delta_alt=alt(i)-alt(first_point) if (delta_alt.ge.delta_alt_limit) then slope_2km=(T(first_point)-T(i))/delta_alt if (slope_2km.lt.slope_limit) then exit_flag=1 else first_point=0 end if end if end if if (first_point.eq.0.and.slope(i).lt.slope_limit) first_point=i end if i=i+i_dir if (i.le.1.or.i.ge.N) exit_flag=1 end do if (first_point.le.0) P_tropo=65.4321 if (first_point.eq.1) P_tropo=64.5432 if (first_point.gt.1) then tmp=(slope_limit-slope(first_point))/(slope(first_point+1)- & & slope(first_point))*(P(first_point+1)-P(first_point)) P_tropo=P(first_point)+tmp ! print*, 'P_tropo= ', tmp, P(first_point), P_tropo end if ! Ajout Marine if (P_tropo .lt. 60. .or. P_tropo .gt. 470.) then P_tropo = 999. endif return end function search_tropopause subroutine cloud_structure(len_cs, rneb_cs, temp_cs, & & emis_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, & & cc_tot_cs, cc_hc_cs, cc_hist_cs, & & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, & & pcld_hc_cs, tcld_hc_cs, & & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, & & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, & & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, & & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, & & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs) INTEGER :: i, n, nss INTEGER, intent(in) :: len_cs REAL, DIMENSION(:), intent(in) :: rneb_cs, temp_cs REAL, DIMENSION(:), intent(in) :: emis_cs, iwco_cs, rad_cs REAL, DIMENSION(:), intent(in) :: pres_cs, dz_cs, rhodz_cs REAL, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, & & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, & & pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, & & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, & & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, & & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, & & em_hist_cs, iwp_hist_cs, & & deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs REAL, DIMENSION(len_cs) :: rneb_ord REAL :: rneb_min REAL, DIMENSION(:), allocatable :: s, s_hc, s_hist, rneb_max REAL, DIMENSION(:), allocatable :: sCb, sThCi, sAnv REAL, DIMENSION(:), allocatable :: iwp_ss, pcld_ss, tcld_ss,& & emis_ss REAL, DIMENSION(:), allocatable :: deltaz_ss, rad_ss CHARACTER (len = 50) :: modname = 'simu_airs.cloud_structure' CHARACTER (len = 160) :: abort_message write(lunout,*) 'dans cloud_structure' call ordonne(len_cs, rneb_cs, rneb_ord) ! Definition des sous_sections rneb_min = rneb_ord(1) nss = 1 do i = 2, size(rneb_cs) if (rneb_ord(i) .gt. rneb_min) then nss = nss+1 rneb_min = rneb_ord(i) endif enddo allocate (s(nss)) allocate (s_hc(nss)) allocate (s_hist(nss)) allocate (rneb_max(nss)) allocate (emis_ss(nss)) allocate (pcld_ss(nss)) allocate (tcld_ss(nss)) allocate (iwp_ss(nss)) allocate (deltaz_ss(nss)) allocate (rad_ss(nss)) allocate (sCb(nss)) allocate (sThCi(nss)) allocate (sAnv(nss)) rneb_min = rneb_ord(1) n = 1 s(1) = rneb_ord(1) s_hc(1) = rneb_ord(1) s_hist(1) = rneb_ord(1) sCb(1) = rneb_ord(1) sThCi(1) = rneb_ord(1) sAnv(1) = rneb_ord(1) rneb_max(1) = rneb_ord(1) do i = 2, size(rneb_cs) if (rneb_ord(i) .gt. rneb_min) then n = n+1 s(n) = rneb_ord(i)-rneb_min s_hc(n) = rneb_ord(i)-rneb_min s_hist(n) = rneb_ord(i)-rneb_min sCb(n) = rneb_ord(i)-rneb_min sThCi(n) = rneb_ord(i)-rneb_min sAnv(n) = rneb_ord(i)-rneb_min rneb_max(n) = rneb_ord(i) rneb_min = rneb_ord(i) endif enddo ! Appel de sous_section do i = 1, nss call sous_section(len_cs, rneb_cs, temp_cs, & & emis_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, & & rneb_max(i),s(i),s_hc(i),s_hist(i), & & sCb(i), sThCi(i), sAnv(i), & & emis_ss(i), & & pcld_ss(i), tcld_ss(i), iwp_ss(i), deltaz_ss(i), rad_ss(i)) enddo ! Caracteristiques de la structure nuageuse cc_tot_cs = 0. cc_hc_cs = 0. cc_hist_cs = 0. cc_Cb_cs = 0. cc_ThCi_cs = 0. cc_Anv_cs = 0. em_hc_cs = 0. iwp_hc_cs = 0. deltaz_hc_cs = 0. em_hist_cs = 0. iwp_hist_cs = 0. deltaz_hist_cs = 0. rad_hist_cs = 0. pcld_hc_cs = 0. tcld_hc_cs = 0. pcld_Cb_cs = 0. tcld_Cb_cs = 0. em_Cb_cs = 0. pcld_ThCi_cs = 0. tcld_ThCi_cs = 0. em_ThCi_cs = 0. pcld_Anv_cs = 0. tcld_Anv_cs = 0. em_Anv_cs = 0. do i = 1, nss cc_tot_cs = cc_tot_cs + s(i) cc_hc_cs = cc_hc_cs + s_hc(i) cc_hist_cs = cc_hist_cs + s_hist(i) cc_Cb_cs = cc_Cb_cs + sCb(i) cc_ThCi_cs = cc_ThCi_cs + sThCi(i) cc_Anv_cs = cc_Anv_cs + sAnv(i) iwp_hc_cs = iwp_hc_cs + s_hc(i)*iwp_ss(i) em_hc_cs = em_hc_cs + s_hc(i)*emis_ss(i) deltaz_hc_cs = deltaz_hc_cs + s_hc(i)*deltaz_ss(i) iwp_hist_cs = iwp_hist_cs + s_hist(i)*iwp_ss(i) em_hist_cs = em_hist_cs + s_hist(i)*emis_ss(i) deltaz_hist_cs = deltaz_hist_cs + s_hist(i)*deltaz_ss(i) rad_hist_cs = rad_hist_cs + s_hist(i)*rad_ss(i) pcld_hc_cs = pcld_hc_cs + s_hc(i)*pcld_ss(i) tcld_hc_cs = tcld_hc_cs + s_hc(i)*tcld_ss(i) pcld_Cb_cs = pcld_Cb_cs + sCb(i)*pcld_ss(i) tcld_Cb_cs = tcld_Cb_cs + sCb(i)*tcld_ss(i) em_Cb_cs = em_Cb_cs + sCb(i)*emis_ss(i) pcld_ThCi_cs = pcld_ThCi_cs + sThCi(i)*pcld_ss(i) tcld_ThCi_cs = tcld_ThCi_cs + sThCi(i)*tcld_ss(i) em_ThCi_cs = em_ThCi_cs + sThCi(i)*emis_ss(i) pcld_Anv_cs = pcld_Anv_cs + sAnv(i)*pcld_ss(i) tcld_Anv_cs = tcld_Anv_cs + sAnv(i)*tcld_ss(i) em_Anv_cs = em_Anv_cs + sAnv(i)*emis_ss(i) enddo deallocate(s) deallocate (s_hc) deallocate (s_hist) deallocate (rneb_max) deallocate (emis_ss) deallocate (pcld_ss) deallocate (tcld_ss) deallocate (iwp_ss) deallocate (deltaz_ss) deallocate (rad_ss) deallocate (sCb) deallocate (sThCi) deallocate (sAnv) call normal_undef(pcld_hc_cs,cc_hc_cs) call normal_undef(tcld_hc_cs,cc_hc_cs) call normal_undef(iwp_hc_cs,cc_hc_cs) call normal_undef(em_hc_cs,cc_hc_cs) call normal_undef(deltaz_hc_cs,cc_hc_cs) call normal_undef(iwp_hist_cs,cc_hist_cs) call normal_undef(em_hist_cs,cc_hist_cs) call normal_undef(deltaz_hist_cs,cc_hist_cs) call normal_undef(rad_hist_cs,cc_hist_cs) call normal_undef(pcld_Cb_cs,cc_Cb_cs) call normal_undef(tcld_Cb_cs,cc_Cb_cs) call normal_undef(em_Cb_cs,cc_Cb_cs) call normal_undef(pcld_ThCi_cs,cc_ThCi_cs) call normal_undef(tcld_ThCi_cs,cc_ThCi_cs) call normal_undef(em_ThCi_cs,cc_ThCi_cs) call normal_undef(pcld_Anv_cs,cc_Anv_cs) call normal_undef(tcld_Anv_cs,cc_Anv_cs) call normal_undef(em_Anv_cs,cc_Anv_cs) ! Tests if (cc_tot_cs .gt. maxval(rneb_cs) .and. & & abs(cc_tot_cs-maxval(rneb_cs)) .gt. 1.e-4 ) then WRITE(abort_message,*) 'cc_tot_cs > max rneb_cs', cc_tot_cs, maxval(rneb_cs) CALL abort_physic(modname,abort_message,1) endif if (iwp_hc_cs .lt. 0.) then abort_message= 'cloud_structure: iwp_hc_cs < 0' CALL abort_physic(modname,abort_message,1) endif end subroutine cloud_structure subroutine normal_undef(num, den) REAL, intent(in) :: den REAL, intent(inout) :: num if (den .ne. 0) then num = num/den else num = undef endif end subroutine normal_undef subroutine normal2_undef(res,num,den) REAL, intent(in) :: den REAL, intent(in) :: num REAL, intent(out) :: res if (den .ne. 0.) then res = num/den else res = undef endif end subroutine normal2_undef subroutine sous_section(len_cs, rneb_cs, temp_cs, & & emis_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, & & rnebmax, stot, shc, shist, & & sCb, sThCi, sAnv, & & emis, pcld, tcld, iwp, deltaz, rad) INTEGER, intent(in) :: len_cs REAL, DIMENSION(len_cs), intent(in) :: rneb_cs, temp_cs REAL, DIMENSION(len_cs), intent(in) :: emis_cs, iwco_cs, & & rneb_ord REAL, DIMENSION(len_cs), intent(in) :: pres_cs, dz_cs, rad_cs REAL, DIMENSION(len_cs), intent(in) :: rhodz_cs REAL, DIMENSION(len_cs) :: tau_cs, w REAL, intent(in) :: rnebmax REAL, intent(inout) :: stot, shc, shist REAL, intent(inout) :: sCb, sThCi, sAnv REAL, intent(out) :: emis, pcld, tcld, iwp, deltaz, rad INTEGER :: i, ideb, ibeg, iend, nuage, visible REAL :: som, som_tau, som_iwc, som_dz, som_rad, tau CHARACTER (len = 50) :: modname = 'simu_airs.sous_section' CHARACTER (len = 160) :: abort_message ! Ponderation: 1 pour les nuages, 0 pour les trous do i = 1, len_cs if (rneb_cs(i) .ge. rnebmax) then w(i) = 1. else w(i) = 0. endif enddo ! Calcul des epaisseurs optiques a partir des emissivites som = 0. do i = 1, len_cs if (emis_cs(i) .eq. 1.) then tau_cs(i) = 10. else tau_cs(i) = -log(1.-emis_cs(i)) endif som = som+tau_cs(i) enddo ideb = 1 nuage = 0 visible = 0 ! Boucle sur les nuages do while (ideb .ne. 0 .and. ideb .le. len_cs) ! Definition d'un nuage de la sous-section call topbot(ideb, w, ibeg, iend) ideb = iend+1 if (ibeg .gt. 0) then nuage = nuage + 1 ! On determine les caracteristiques du nuage ! (ep. optique, ice water path, pression, temperature) call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, & & som_tau, som_iwc, som_dz, som_rad) ! On masque le nuage s'il n'est pas detectable call masque(ibeg, iend, som_tau, visible, w) endif ! Fin boucle nuages enddo ! Analyse du nuage detecte call topbot(1, w, ibeg, iend) if (ibeg .gt. 0) then call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, & & som_tau, som_iwc, som_dz, som_rad) tau = som_tau emis = 1. - exp(-tau) iwp = som_iwc deltaz = som_dz rad = som_rad if (pcld .gt. p_thresh) then shc = 0. shist = 0. sCb = 0. sThCi = 0. sAnv = 0. else if (emis .lt. em_min .or. emis .gt. em_max & & .or. tcld .gt. 230.) then shist = 0. endif if (emis .lt. 0.98) then sCb = 0. endif if (emis .gt. 0.5 .or. emis .lt. 0.1) then sThCi = 0. endif if (emis .le. 0.5 .or. emis .ge. 0.98) then sAnv = 0. endif endif else tau = 0. emis = 0. iwp = 0. deltaz = 0. pcld = 0. tcld = 0. stot = 0. shc = 0. shist = 0. rad = 0. sCb = 0. sThCi = 0. sAnv = 0. endif ! Tests if (iwp .lt. 0.) then WRITE(abort_message,*) 'ideb iwp =', ideb, iwp CALL abort_physic(modname,abort_message,1) endif if (deltaz .lt. 0.) then WRITE(abort_message,*)'ideb deltaz =', ideb, deltaz CALL abort_physic(modname,abort_message,1) endif if (emis .lt. 0.048 .and. emis .ne. 0.) then WRITE(abort_message,*) 'ideb emis =', ideb, emis CALL abort_physic(modname,abort_message,1) endif end subroutine sous_section subroutine masque (ibeg, iend, som_tau, & & visible, w) INTEGER, intent(in) :: ibeg, iend REAL, intent(in) :: som_tau INTEGER, intent(inout) :: visible REAL, DIMENSION(:), intent(inout) :: w INTEGER :: i ! Masque ! Cas ou il n'y a pas de nuage visible au-dessus if (visible .eq. 0) then if (som_tau .lt. tau_thresh) then do i = ibeg, iend w(i) = 0. enddo else visible = 1 endif ! Cas ou il y a un nuage visible au-dessus else do i = ibeg, iend w(i) = 0. enddo endif end subroutine masque subroutine caract (ibeg, iend, temp_cs, tau_cs, iwco_cs, & & pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, & & som_tau, som_iwc, som_dz, som_rad) INTEGER, intent(in) :: ibeg, iend REAL, DIMENSION(:), intent(in) :: tau_cs, iwco_cs, temp_cs REAL, DIMENSION(:), intent(in) :: pres_cs, dz_cs, rad_cs REAL, DIMENSION(:), intent(in) :: rhodz_cs REAL, intent(out) :: som_tau, som_iwc, som_dz, som_rad REAL , intent(out) :: pcld, tcld INTEGER :: i, ibase, imid CHARACTER (len = 50) :: modname = 'simu_airs.caract' CHARACTER (len = 160) :: abort_message ! Somme des epaisseurs optiques et des contenus en glace sur le nuage som_tau = 0. som_iwc = 0. som_dz = 0. som_rad = 0. ibase = -100 do i = ibeg, iend som_tau = som_tau + tau_cs(i) som_dz = som_dz + dz_cs(i) som_iwc = som_iwc + iwco_cs(i)*1000*rhodz_cs(i) ! en g/m2 som_rad = som_rad + rad_cs(i)*dz_cs(i) if (som_tau .gt. 3. .and. ibase .eq. -100) then ibase = i-1 endif enddo if (som_dz .ne. 0.) then som_rad = som_rad/som_dz else write(*,*) 'som_dez = 0 STOP' write(*,*) 'ibeg, iend =', ibeg, iend do i = ibeg, iend write(*,*) dz_cs(i), rhodz_cs(i) enddo abort_message='see above' CALL abort_physic(modname,abort_message,1) endif ! Determination de Pcld et Tcld if (ibase .lt. ibeg) then ibase = ibeg endif imid = (ibeg+ibase)/2 pcld = pres_cs(imid)/100. ! pcld en hPa tcld = temp_cs(imid) end subroutine caract subroutine topbot(ideb,w,ibeg,iend) INTEGER, intent(in) :: ideb REAL, DIMENSION(:), intent(in) :: w INTEGER, intent(out) :: ibeg, iend INTEGER :: i, itest itest = 0 ibeg = 0 iend = 0 do i = ideb, size(w) if (w(i) .eq. 1. .and. itest .eq. 0) then ibeg = i itest = 1 endif enddo i = ibeg do while (w(i) .eq. 1. .and. i .le. size(w)) i = i+1 enddo iend = i-1 end subroutine topbot subroutine ordonne(len_cs, rneb_cs, rneb_ord) INTEGER, intent(in) :: len_cs REAL, DIMENSION(:), intent(in) :: rneb_cs REAL, DIMENSION(:), intent(out) :: rneb_ord INTEGER :: i, j, ind_min REAL, DIMENSION(len_cs) :: rneb REAL :: rneb_max do i = 1, size(rneb_cs) rneb(i) = rneb_cs(i) enddo do j = 1, size(rneb_cs) rneb_max = 100. do i = 1, size(rneb_cs) if (rneb(i) .lt. rneb_max) then rneb_max = rneb(i) ind_min = i endif enddo rneb_ord(j) = rneb_max rneb(ind_min) = 100. enddo end subroutine ordonne subroutine sim_mesh(rneb_1D, temp_1D, emis_1D, & & iwcon_1D, rad_1D, & & pres, dz, & & rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh,& & tcld_hc_mesh, & & em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, & & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, & & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, & & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, & & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, & & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh) USE dimphy REAL, DIMENSION(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, & & iwcon_1D, rad_1D REAL, DIMENSION(klev), intent(in) :: pres, dz, rhodz_1D REAL, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh REAL, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh REAL, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, & & iwp_hc_mesh REAL, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh REAL, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, & & em_ThCi_mesh REAL, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh REAL, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh REAL, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh REAL, DIMENSION(:), allocatable :: rneb_cs, temp_cs, emis_cs, & & iwco_cs REAL, DIMENSION(:), allocatable :: pres_cs, dz_cs, rad_cs, & & rhodz_cs INTEGER :: i,j,l INTEGER :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics REAL :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,& & som_hist REAL :: som_emi_hist, som_iwp_hist, som_deltaz_hc, & & som_deltaz_hist REAL :: som_rad_hist REAL :: som_Cb, som_ThCi, som_Anv REAL :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb REAL :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv REAL :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi REAL :: tsom_tot, tsom_hc, tsom_hist REAL :: prod, prod_hh REAL :: cc_tot_cs, cc_hc_cs, cc_hist_cs REAL :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs REAL :: pcld_hc_cs, tcld_hc_cs REAL :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs REAL :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs REAL :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs REAL :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs REAL :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs REAL, DIMENSION(klev) :: test_tot, test_hc, test_hist REAL, DIMENSION(klev) :: test_pcld, test_tcld, test_em, test_iwp CHARACTER (len = 50) :: modname = 'simu_airs.sim_mesh' CHARACTER (len = 160) :: abort_message do j = 1, klev WRITE(lunout,*) 'simu_airs, j, rneb_1D =', rneb_1D(j) enddo ! Definition des structures nuageuses, de la plus haute a la plus basse num_cs = 0 ltop = klev-1 prod = 1. som_emi_hc = 0. som_emi_hist = 0. som_pcl_hc = 0. som_tcl_hc = 0. som_iwp_hc = 0. som_iwp_hist = 0. som_deltaz_hc = 0. som_deltaz_hist = 0. som_rad_hist = 0. som_hc = 0. som_hist = 0. som_Cb = 0. som_ThCi = 0. som_Anv = 0. som_emi_Cb = 0. som_pcld_Cb = 0. som_tcld_Cb = 0. som_emi_ThCi = 0. som_pcld_ThCi = 0. som_tcld_ThCi = 0. som_emi_Anv = 0. som_pcld_Anv = 0. som_tcld_Anv = 0. tsom_tot = 0. tsom_hc = 0. tsom_hist = 0. do while (ltop .ge. 1) ! Boucle sur toute la colonne itop = 0 do l = ltop,1,-1 if (itop .eq. 0 .and. rneb_1D(l) .gt. 0.001 ) then itop = l endif enddo ibot = itop do while (rneb_1D(ibot) .gt. 0.001 .and. ibot .ge. 1) ibot = ibot -1 enddo ibot = ibot+1 if (itop .gt. 0) then ! itop > 0 num_cs = num_cs +1 len_cs = itop-ibot+1 ! Allocation et definition des variables de la structure nuageuse ! le premier indice denote ce qui est le plus haut allocate (rneb_cs(len_cs)) allocate (temp_cs(len_cs)) allocate (emis_cs(len_cs)) allocate (iwco_cs(len_cs)) allocate (pres_cs(len_cs)) allocate (dz_cs(len_cs)) allocate (rad_cs(len_cs)) allocate (rhodz_cs(len_cs)) ics = 0 do i = itop, ibot, -1 ics = ics + 1 rneb_cs(ics) = rneb_1D(i) temp_cs(ics) = temp_1D(i) emis_cs(ics) = emis_1D(i) iwco_cs(ics) = iwcon_1D(i) rad_cs(ics) = rad_1D(i) pres_cs(ics) = pres(i) dz_cs(ics) = dz(i) rhodz_cs(ics) = rhodz_1D(i) enddo ! Appel du sous_programme cloud_structure call cloud_structure(len_cs,rneb_cs,temp_cs,emis_cs,iwco_cs,& & pres_cs, dz_cs, rhodz_cs, rad_cs, & & cc_tot_cs, cc_hc_cs, cc_hist_cs, & & cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, & & pcld_hc_cs, tcld_hc_cs, & & em_hc_cs, iwp_hc_cs, deltaz_hc_cs, & & pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, & & pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, & & pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, & & em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs) deallocate (rneb_cs) deallocate (temp_cs) deallocate (emis_cs) deallocate (iwco_cs) deallocate (pres_cs) deallocate (dz_cs) deallocate (rad_cs) deallocate (rhodz_cs) ! Pour la couverture nuageuse sur la maille prod_hh = prod prod = prod*(1.-cc_tot_cs) ! Pour les autres variables definies sur la maille som_emi_hc = som_emi_hc + em_hc_cs*cc_hc_cs*prod_hh som_iwp_hc = som_iwp_hc + iwp_hc_cs*cc_hc_cs*prod_hh som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs*cc_hc_cs*prod_hh som_emi_Cb = som_emi_Cb + em_Cb_cs*cc_Cb_cs*prod_hh som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs*cc_Cb_cs*prod_hh som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs*cc_Cb_cs*prod_hh som_emi_ThCi = som_emi_ThCi + em_ThCi_cs*cc_ThCi_cs*prod_hh som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs*cc_ThCi_cs*prod_hh som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs*cc_ThCi_cs*prod_hh som_emi_Anv = som_emi_Anv + em_Anv_cs*cc_Anv_cs*prod_hh som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs*cc_Anv_cs*prod_hh som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs*cc_Anv_cs*prod_hh som_emi_hist = som_emi_hist + em_hist_cs*cc_hist_cs*prod_hh som_iwp_hist = som_iwp_hist + iwp_hist_cs*cc_hist_cs*prod_hh som_deltaz_hist = som_deltaz_hist + & & deltaz_hist_cs*cc_hist_cs*prod_hh som_rad_hist = som_rad_hist + rad_hist_cs*cc_hist_cs*prod_hh som_pcl_hc = som_pcl_hc + pcld_hc_cs*cc_hc_cs*prod_hh som_tcl_hc = som_tcl_hc + tcld_hc_cs*cc_hc_cs*prod_hh som_hc = som_hc + cc_hc_cs*prod_hh som_hist = som_hist + cc_hist_cs*prod_hh som_Cb = som_Cb + cc_Cb_cs*prod_hh som_ThCi = som_ThCi + cc_ThCi_cs*prod_hh som_Anv = som_Anv + cc_Anv_cs*prod_hh ! Pour test call test_bornes('cc_tot_cs ',cc_tot_cs,1.,0.) call test_bornes('cc_hc_cs ',cc_hc_cs,1.,0.) call test_bornes('cc_hist_cs ',cc_hist_cs,1.,0.) call test_bornes('pcld_hc_cs ',pcld_hc_cs,1200.,0.) call test_bornes('tcld_hc_cs ',tcld_hc_cs,1000.,100.) call test_bornes('em_hc_cs ',em_hc_cs,1000.,0.048) test_tot(num_cs) = cc_tot_cs test_hc(num_cs) = cc_hc_cs test_hist(num_cs) = cc_hist_cs test_pcld(num_cs) = pcld_hc_cs test_tcld(num_cs) = tcld_hc_cs test_em(num_cs) = em_hc_cs test_iwp(num_cs) = iwp_hc_cs tsom_tot = tsom_tot + cc_tot_cs tsom_hc = tsom_hc + cc_hc_cs tsom_hist = tsom_hist + cc_hist_cs endif ! itop > 0 ltop = ibot -1 enddo ! fin de la boucle sur la colonne N_CS = num_cs ! Determination des variables de sortie if (N_CS .gt. 0) then ! if N_CS>0 cc_tot_mesh = 1. - prod cc_hc_mesh = som_hc cc_hist_mesh = som_hist cc_Cb_mesh = som_Cb cc_ThCi_mesh = som_ThCi cc_Anv_mesh = som_Anv call normal2_undef(pcld_hc_mesh,som_pcl_hc, & & cc_hc_mesh) call normal2_undef(tcld_hc_mesh,som_tcl_hc, & & cc_hc_mesh) call normal2_undef(em_hc_mesh,som_emi_hc, & & cc_hc_mesh) call normal2_undef(iwp_hc_mesh,som_iwp_hc, & & cc_hc_mesh) call normal2_undef(deltaz_hc_mesh,som_deltaz_hc, & & cc_hc_mesh) call normal2_undef(em_Cb_mesh,som_emi_Cb, & & cc_Cb_mesh) call normal2_undef(tcld_Cb_mesh,som_tcld_Cb, & & cc_Cb_mesh) call normal2_undef(pcld_Cb_mesh,som_pcld_Cb, & & cc_Cb_mesh) call normal2_undef(em_ThCi_mesh,som_emi_ThCi, & & cc_ThCi_mesh) call normal2_undef(tcld_ThCi_mesh,som_tcld_ThCi, & & cc_ThCi_mesh) call normal2_undef(pcld_ThCi_mesh,som_pcld_ThCi, & & cc_ThCi_mesh) call normal2_undef(em_Anv_mesh,som_emi_Anv, & & cc_Anv_mesh) call normal2_undef(tcld_Anv_mesh,som_tcld_Anv, & & cc_Anv_mesh) call normal2_undef(pcld_Anv_mesh,som_pcld_Anv, & & cc_Anv_mesh) call normal2_undef(em_hist_mesh,som_emi_hist, & & cc_hist_mesh) call normal2_undef(iwp_hist_mesh,som_iwp_hist, & & cc_hist_mesh) call normal2_undef(deltaz_hist_mesh,som_deltaz_hist, & & cc_hist_mesh) call normal2_undef(rad_hist_mesh,som_rad_hist, & & cc_hist_mesh) ! Tests ! Tests if (cc_tot_mesh .gt. tsom_tot .and. & & abs(cc_tot_mesh-tsom_tot) .gt. 1.e-4) then WRITE(abort_message,*)'cc_tot_mesh > tsom_tot', cc_tot_mesh, tsom_tot CALL abort_physic(modname,abort_message,1) endif if (cc_tot_mesh .lt. maxval(test_tot(1:N_CS)) .and. & & abs(cc_tot_mesh-maxval(test_tot(1:N_CS))) .gt. 1.e-4) then WRITE(abort_message,*) 'cc_tot_mesh < max', cc_tot_mesh, maxval(test_tot(1:N_CS)) CALL abort_physic(modname,abort_message,1) endif if (cc_hc_mesh .gt. tsom_hc .and. & & abs(cc_hc_mesh-tsom_hc) .gt. 1.e-4) then WRITE(abort_message,*) 'cc_hc_mesh > tsom_hc', cc_hc_mesh, tsom_hc CALL abort_physic(modname,abort_message,1) endif if (cc_hc_mesh .lt. maxval(test_hc(1:N_CS)) .and. & & abs(cc_hc_mesh-maxval(test_hc(1:N_CS))) .gt. 1.e-4) then WRITE(abort_message,*) 'cc_hc_mesh < max', cc_hc_mesh, maxval(test_hc(1:N_CS)) CALL abort_physic(modname,abort_message,1) endif if (cc_hist_mesh .gt. tsom_hist .and. & & abs(cc_hist_mesh-tsom_hist) .gt. 1.e-4) then WRITE(abort_message,*) 'cc_hist_mesh > tsom_hist', cc_hist_mesh, tsom_hist CALL abort_physic(modname,abort_message,1) endif if (cc_hist_mesh .lt. 0.) then WRITE(abort_message,*) 'cc_hist_mesh < 0', cc_hist_mesh CALL abort_physic(modname,abort_message,1) endif if ((pcld_hc_mesh .gt. maxval(test_pcld(1:N_CS)) .or. & & pcld_hc_mesh .lt. minval(test_pcld(1:N_CS))) .and. & & abs(pcld_hc_mesh-maxval(test_pcld(1:N_CS))) .gt. 1. .and. & & maxval(test_pcld(1:N_CS)) .ne. 999. & & .and. minval(test_pcld(1:N_CS)) .ne. 999.) then WRITE(abort_message,*) 'pcld_hc_mesh est faux', pcld_hc_mesh, maxval(test_pcld(1:N_CS)), & & minval(test_pcld(1:N_CS)) CALL abort_physic(modname,abort_message,1) endif if ((tcld_hc_mesh .gt. maxval(test_tcld(1:N_CS)) .or. & & tcld_hc_mesh .lt. minval(test_tcld(1:N_CS))) .and. & & abs(tcld_hc_mesh-maxval(test_tcld(1:N_CS))) .gt. 0.1 .and. & & maxval(test_tcld(1:N_CS)) .ne. 999. & & .and. minval(test_tcld(1:N_CS)) .ne. 999.) then WRITE(abort_message,*) 'tcld_hc_mesh est faux', tcld_hc_mesh, maxval(test_tcld(1:N_CS)), & & minval(test_tcld(1:N_CS)) CALL abort_physic(modname,abort_message,1) endif if ((em_hc_mesh .gt. maxval(test_em(1:N_CS)) .or. & & em_hc_mesh .lt. minval(test_em(1:N_CS))) .and. & & abs(em_hc_mesh-maxval(test_em(1:N_CS))) .gt. 1.e-4 .and. & & minval(test_em(1:N_CS)) .ne. 999. .and. & & maxval(test_em(1:N_CS)) .ne. 999. ) then WRITE(abort_message,*) 'em_hc_mesh est faux', em_hc_mesh, maxval(test_em(1:N_CS)), & & minval(test_em(1:N_CS)) CALL abort_physic(modname,abort_message,1) endif else ! if N_CS>0 cc_tot_mesh = 0. cc_hc_mesh = 0. cc_hist_mesh = 0. cc_Cb_mesh = 0. cc_ThCi_mesh = 0. cc_Anv_mesh = 0. iwp_hc_mesh = undef deltaz_hc_mesh = undef em_hc_mesh = undef iwp_hist_mesh = undef deltaz_hist_mesh = undef rad_hist_mesh = undef em_hist_mesh = undef pcld_hc_mesh = undef tcld_hc_mesh = undef pcld_Cb_mesh = undef tcld_Cb_mesh = undef em_Cb_mesh = undef pcld_ThCi_mesh = undef tcld_ThCi_mesh = undef em_ThCi_mesh = undef pcld_Anv_mesh = undef tcld_Anv_mesh = undef em_Anv_mesh = undef endif ! if N_CS>0 end subroutine sim_mesh subroutine test_bornes(sx,x,bsup,binf) REAL, intent(in) :: x, bsup, binf character*14, intent(in) :: sx CHARACTER (len = 50) :: modname = 'simu_airs.test_bornes' CHARACTER (len = 160) :: abort_message if (x .gt. bsup .or. x .lt. binf) then WRITE(abort_message,*) sx, 'est faux', sx, x CALL abort_physic(modname,abort_message,1) endif end subroutine test_bornes end module m_simu_airs subroutine simu_airs & & (itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, & & geop_airs, pplay_airs, paprs_airs, & & map_prop_hc,map_prop_hist,& & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,& & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,& & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,& & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,& & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,& & map_ntot,map_hc,map_hist,& & map_Cb,map_ThCi,map_Anv,alt_tropo ) USE dimphy USE m_simu_airs USE yomcst_mod_h IMPLICIT NONE INTEGER,intent(in) :: itap REAL, DIMENSION(klon,klev), intent(in) :: & & rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, & & rad_airs, geop_airs, pplay_airs, paprs_airs REAL, DIMENSION(klon,klev) :: & & rhodz_airs, rho_airs, iwcon_airs REAL, DIMENSION(klon),intent(out) :: alt_tropo REAL, DIMENSION(klev) :: rneb_1D, temp_1D, & & emis_1D, rad_1D, pres_1D, alt_1D, & & rhodz_1D, dz_1D, iwcon_1D INTEGER :: i, j REAL :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh REAL :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh REAL :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh REAL :: em_hist_mesh, iwp_hist_mesh REAL :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh REAL :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh REAL :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh REAL :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh REAL, DIMENSION(klon),intent(out) :: map_prop_hc, map_prop_hist REAL, DIMENSION(klon),intent(out) :: map_emis_hc, map_iwp_hc REAL, DIMENSION(klon),intent(out) :: map_deltaz_hc, map_pcld_hc REAL, DIMENSION(klon),intent(out) :: map_tcld_hc REAL, DIMENSION(klon),intent(out) :: map_emis_Cb,map_pcld_Cb,map_tcld_Cb REAL, DIMENSION(klon),intent(out) :: & & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi REAL, DIMENSION(klon),intent(out) :: & & map_emis_Anv,map_pcld_Anv,map_tcld_Anv REAL, DIMENSION(klon),intent(out) :: & & map_emis_hist,map_iwp_hist,map_deltaz_hist,& & map_rad_hist REAL, DIMENSION(klon),intent(out) :: map_ntot,map_hc,map_hist REAL, DIMENSION(klon),intent(out) :: map_Cb,map_ThCi,map_Anv write(*,*) 'simu_airs' write(*,*) 'itap, klon, klev', itap, klon, klev write(*,*) 'RG, RD =', RG, RD ! Definition des variables 1D do i = 1, klon do j = 1, klev-1 rhodz_airs(i,j) = & & (paprs_airs(i,j)-paprs_airs(i,j+1))/RG enddo rhodz_airs(i,klev) = 0. enddo do i = 1, klon do j = 1,klev rho_airs(i,j) = & & pplay_airs(i,j)/(temp_airs(i,j)*RD) if (rneb_airs(i,j) .gt. 0.001) then iwcon_airs(i,j) = iwcon0_airs(i,j)/rneb_airs(i,j) else iwcon_airs(i,j) = 0. endif enddo enddo !============================================================================= do i = 1, klon ! boucle sur les points de grille !============================================================================= do j = 1,klev rneb_1D(j) = rneb_airs(i,j) temp_1D(j) = temp_airs(i,j) emis_1D(j) = cldemi_airs(i,j) iwcon_1D(j) = iwcon_airs(i,j) rad_1D(j) = rad_airs(i,j) pres_1D(j) = pplay_airs(i,j) alt_1D(j) = geop_airs(i,j)/RG rhodz_1D(j) = rhodz_airs(i,j) dz_1D(j) = rhodz_airs(i,j)/rho_airs(i,j) enddo alt_tropo(i) = & & search_tropopause(pres_1D/100.,temp_1D,alt_1D,klev) ! Appel du ss-programme sim_mesh ! if (itap .eq. 1 ) then call sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, & & pres_1D, dz_1D, rhodz_1D, & & cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, & & pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, & & deltaz_hc_mesh,& & cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, & & pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, & & pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, & & pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, & & em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh) write(*,*) '====================================' write(*,*) 'itap, i:', itap, i write(*,*) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, & & iwp_hc, em_hist, iwp_hist =' write(*,*) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh write(*,*) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh write(*,*) em_hist_mesh, iwp_hist_mesh ! endif ! Definition des variables a ecrire dans le fichier de sortie call normal2_undef(map_prop_hc(i),cc_hc_mesh, & & cc_tot_mesh) call normal2_undef(map_prop_hist(i),cc_hist_mesh, & & cc_tot_mesh) map_emis_hc(i) = em_hc_mesh map_iwp_hc(i) = iwp_hc_mesh map_deltaz_hc(i) = deltaz_hc_mesh map_pcld_hc(i) = pcld_hc_mesh map_tcld_hc(i) = tcld_hc_mesh map_emis_Cb(i) = em_Cb_mesh map_pcld_Cb(i) = pcld_Cb_mesh map_tcld_Cb(i) = tcld_Cb_mesh map_emis_ThCi(i) = em_ThCi_mesh map_pcld_ThCi(i) = pcld_ThCi_mesh map_tcld_ThCi(i) = tcld_ThCi_mesh map_emis_Anv(i) = em_Anv_mesh map_pcld_Anv(i) = pcld_Anv_mesh map_tcld_Anv(i) = tcld_Anv_mesh map_emis_hist(i) = em_hist_mesh map_iwp_hist(i) = iwp_hist_mesh map_deltaz_hist(i) = deltaz_hist_mesh map_rad_hist(i) = rad_hist_mesh map_ntot(i) = cc_tot_mesh map_hc(i) = cc_hc_mesh map_hist(i) = cc_hist_mesh map_Cb(i) = cc_Cb_mesh map_ThCi(i) = cc_ThCi_mesh map_Anv(i) = cc_Anv_mesh enddo ! fin boucle sur les points de grille end subroutine simu_airs