Changeset 2669 for LMDZ5/branches/testing
- Timestamp:
- Oct 14, 2016, 2:57:28 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 3 deleted
- 33 edited
- 6 copied
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2642-2664
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/DefLists/context_lmdz.xml
r2594 r2669 31 31 <grid_group id="vertical" axis_ref="presnivs" /> 32 32 </grid_> --> 33 33 34 <!-- Define Scalar grid for GHG, orbital parameters and solar constants --> 35 <grid_definition> 36 <grid id="grid_scalar" /> 37 </grid_definition> 38 34 39 <!-- Define groups of vertical axes --> 35 40 <axis_definition> -
LMDZ5/branches/testing/DefLists/field_def_lmdz.xml
r2641 r2669 867 867 </field_group> 868 868 869 870 <field_group id="field_scalar" operation="instant" freq_op="30d" grid_ref="grid_scalar"> 871 <field id="R_ecc" long_name="R_ecc" unit="-" /> 872 <field id="R_peri" long_name="R_peri" unit="-" /> 873 <field id="R_incl" long_name="R_incl" unit="deg" /> 874 <field id="solaire" long_name="solaire" unit="W/m2" /> 875 876 <field id="rsun1" long_name="Fraction constante solaire bande 1" unit="W/m2" /> 877 <field id="rsun2" long_name="Fraction constante solaire bande 2" unit="W/m2" /> 878 <field id="rsun3" long_name="Fraction constante solaire bande 3" unit="W/m2" /> 879 <field id="rsun4" long_name="Fraction constante solaire bande 4" unit="W/m2" /> 880 <field id="rsun5" long_name="Fraction constante solaire bande 5" unit="W/m2" /> 881 <field id="rsun6" long_name="Fraction constante solaire bande 6" unit="W/m2" /> 882 883 <field id="co2_ppm" long_name="co2_ppm" unit="ppm" /> 884 <field id="CH4_ppb" long_name="CH4_ppb" unit="ppb" /> 885 <field id="N2O_ppb" long_name="N2O_ppb" unit="ppb" /> 886 <field id="CFC11_ppt" long_name="CFC11_ppt" unit="ppt" /> 887 <field id="CFC12_ppt" long_name="CFC12_ppt" unit="ppt" /> 888 </field_group> 889 869 890 </field_definition> 870 -
LMDZ5/branches/testing/DefLists/file_def_histmth_lmdz.xml
r2594 r2669 10 10 <!-- </field_group> --> 11 11 12 13 <field_group operation="instant" freq_op="30d"> 14 <field field_ref="R_ecc" level="1" name="R_ecc" /> 15 <field field_ref="R_peri" level="1" name="R_peri" /> 16 <field field_ref="R_incl" level="1" name="R_incl" /> 17 <field field_ref="solaire" level="1" name="solaire" /> 18 19 <field field_ref="rsun1" level="1" name="rsun1" /> 20 <field field_ref="rsun2" level="1" name="rsun2" /> 21 <field field_ref="rsun3" level="1" name="rsun3" /> 22 <field field_ref="rsun4" level="1" name="rsun4" /> 23 <field field_ref="rsun5" level="1" name="rsun5" /> 24 <field field_ref="rsun6" level="1" name="rsun6" /> 25 26 <field field_ref="co2_ppm" level="1" name="co2_ppm" /> 27 <field field_ref="CH4_ppb" level="1" name="CH4_ppb" /> 28 <field field_ref="N2O_ppb" level="1" name="N2O_ppb" /> 29 <field field_ref="CFC11_ppt" level="1" name="CFC11_ppt" /> 30 <field field_ref="CFC12_ppt" level="1" name="CFC12_ppt" /> 31 32 </field_group> 33 12 34 <!-- VARS 2D --> 13 35 <field_group operation="average" freq_op="1ts"> -
LMDZ5/branches/testing/DefLists/gcm.def_144x142
r2471 r2669 1 ## $Id : gcm.def 1403 2010-07-01 09:02:53Z fairhead$1 ## $Id$ 2 2 ## nombre de pas par jour (multiple de iperiod) ( ici pour dt = 1 min ) 3 day_step= 9603 day_step=720 4 4 ## periode pour le pas Matsuno (en pas) 5 5 iperiod=5 … … 13 13 niterh=2 14 14 ## temps de dissipation des plus petites long.d ondes pour u,v (gradiv) 15 tetagdiv= 1800.15 tetagdiv=3600. 16 16 ## temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot) 17 #tetagrot=5400.18 17 tetagrot=7200. 19 18 ## temps de dissipation des plus petites long.d ondes pour h ( divgrad) 20 tetatemp= 5400.19 tetatemp=7200. 21 20 ## coefficient pour gamdissip 22 21 coefdis=0. … … 33 32 read_start=y 34 33 ## periode de la physique (en pas dynamiques, n'a de sens que si iflag_phys=1) 34 ## iphysiq est modifie pour la nouvelle physique par le lmdz.driver pour les autres cas le valeur par defaut est utilise 35 35 iphysiq=_AUTO_ : DEFAULT = 5 36 #37 #38 ### Parametres du zoom39 36 ## longitude en degres du centre du zoom 40 37 clon=0. … … 57 54 ## Fonction f(y) avec y = Sin(latit.) si = .true. , sinon y = latit. 58 55 ysinus=y 56 ## 57 ngroup=4 -
LMDZ5/branches/testing/DefLists/iodef_dev.xml
r2056 r2669 1 1 <?xml version="1.0"?> 2 2 <simulation> 3 <context id="LMDZ" calendar_type="D360" start_date="1980-01-01 00:00:00"> 3 <context id="LMDZ"> 4 5 <!-- <calendar type="D360" start_date="1980-01-01 00:00:00" /> --> 6 4 7 <!-- Definition of model variables --> 5 8 <field_definition prec="4" … … 16 19 <!-- 3D variables --> 17 20 <field_group id="fields_3D" 18 domain_ref="dom_glo" >21 domain_ref="dom_glo" axis_ref="presnivs"> 19 22 <field id="temperature" 20 23 long_name="Atmospheric temperature" … … 39 42 name="Xhistins" 40 43 output_freq="12ts" 44 type="one_file" 41 45 enabled=".TRUE."> <!-- NB: output_freq in physics ts--> 42 46 … … 46 50 <field field_ref="ps" /> 47 51 </field_group> 48 <!-- <field_group field_group_ref="fields_2D"49 operation="instant" /> doesn't work ?!? -->50 52 51 53 <!-- VARS 3D --> … … 58 60 <field field_ref="v" /> 59 61 </field_group> 60 <!-- <field_group field_group_ref="fields_3D"61 operation="instant" /> doesn't work ?!? -->62 62 </file> 63 63 … … 65 65 name="Xdiurnalave" 66 66 output_freq="1d" 67 type="one_file" 67 68 enabled=".true."> 68 69 <!-- VARS 2D --> … … 93 94 wxios_add_vaxis --> 94 95 <axis_definition> 95 <!--96 <axis_group id="presnivs"97 standard_name="Pseudo-pressure of model vertical levels"98 unit="Pa">99 </axis_group>100 -->101 96 <axis id="presnivs" 102 97 standard_name="Pseudo-pressure of model vertical levels" … … 114 109 115 110 <variable_group id="parameters" > 116 <variable id="using_server" type="bool ean">true</variable>111 <variable id="using_server" type="bool">false</variable> 117 112 <variable id="info_level" type="int">10</variable> 118 113 </variable_group> -
LMDZ5/branches/testing/DefLists/physiq.def_NPv5.65
r2594 r2669 36 36 gwd_rando_sat=0.2 37 37 38 ok_strato=y39 ok_hines=n40 38 sso_gkdrag=0.2 41 39 sso_gfrcri=0.70 -
LMDZ5/branches/testing/DefLists/physiq.def_NPv5.67
r2594 r2669 36 36 gwd_rando_sat=0.2 37 37 38 ok_strato=y39 ok_hines=n40 38 sso_gkdrag=0.2 41 39 sso_gfrcri=0.70 -
LMDZ5/branches/testing/libf/dynphy_lonlat/phydev/callphysiq_mod.F90
r2435 r2669 12 12 jD_cur,jH_cur_split,zdt_split, & 13 13 zplev_omp,zplay_omp, & 14 zp hi_omp,zphis_omp,&14 zpk_omp,zphi_omp,zphis_omp, & 15 15 presnivs_omp, & 16 16 zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp, & … … 34 34 REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa) 35 35 REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa) 36 REAL,INTENT(IN) :: zpk_omp(klon,llm) ! Exner function 36 37 REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer 37 38 REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential -
LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
r2641 r2669 111 111 INTEGER :: flag_aerosol 112 112 INTEGER :: flag_aerosol_strat 113 LOGICAL :: flag_bc_internal_mixture 113 114 LOGICAL :: new_aod 114 115 REAL :: bl95_b0, bl95_b1 … … 131 132 ok_ade, ok_aie, ok_cdnc, aerosol_couple, & 132 133 flag_aerosol, flag_aerosol_strat, new_aod, & 133 bl95_b0, bl95_b1,&134 flag_bc_internal_mixture, bl95_b0, bl95_b1, & 134 135 read_climoz, & 135 136 alp_offset) -
LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
r2641 r2669 9 9 nbp, communicator, & 10 10 punjours, pdayref,ptimestep, & 11 rlatu ,rlatv,rlonu,rlonv,aire,cu,cv,&11 rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn, & 12 12 prad,pg,pr,pcpp,iflag_phys) 13 13 USE dimphy, ONLY: init_dimphy … … 44 44 USE mod_phys_lmdz_omp_data, ONLY: klon_omp 45 45 #endif 46 USE ioipsl_getin_p_mod, ONLY: getin_p 47 USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom 46 48 IMPLICIT NONE 47 49 … … 52 54 53 55 include "dimensions.h" 56 include "paramet.h" 54 57 include "iniprint.h" 55 58 include "tracstoke.h" 59 include "comgeom.h" 56 60 57 61 REAL, INTENT (IN) :: prad ! radius of the planet (m) … … 65 69 INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process 66 70 INTEGER, INTENT(IN) :: communicator ! MPI communicator 67 REAL, INTENT (IN) :: rlatu (jj+1) ! latitudes of the physics grid68 REAL, INTENT (IN) :: rlatv (jj) ! latitude boundaries of the physics grid69 REAL, INTENT (IN) :: rlonv (ii+1) ! longitudes of the physics grid70 REAL, INTENT (IN) :: rlonu (ii+1) ! longitude boundaries of the physics grid71 REAL, INTENT (IN) :: aire (ii+1,jj+1) ! area of the dynamics grid (m2)72 REAL, INTENT (IN) :: cu ((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)73 REAL, INTENT (IN) :: cv ((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)71 REAL, INTENT (IN) :: rlatudyn(jj+1) ! latitudes of the physics grid 72 REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid 73 REAL, INTENT (IN) :: rlonvdyn(ii+1) ! longitudes of the physics grid 74 REAL, INTENT (IN) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid 75 REAL, INTENT (IN) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2) 76 REAL, INTENT (IN) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u) 77 REAL, INTENT (IN) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v) 74 78 INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation 75 79 REAL, INTENT (IN) :: ptimestep !physics time step (s) … … 80 84 CHARACTER (LEN=20) :: modname = 'iniphysiq' 81 85 CHARACTER (LEN=80) :: abort_message 82 86 87 LOGICAL :: slab_hdiff 88 INTEGER :: slab_ekman 89 CHARACTER (LEN = 6) :: type_ocean 83 90 84 91 #ifndef CPP_PARA … … 92 99 CALL inigeomphy(ii,jj,nlayer, & 93 100 nbp, communicator, & 94 rlatu ,rlatv, &95 rlonu ,rlonv, &96 aire ,cu,cv)101 rlatudyn,rlatvdyn, & 102 rlonudyn,rlonvdyn, & 103 airedyn,cudyn,cvdyn) 97 104 98 105 ! --> now initialize things specific to the phylmd physics package … … 117 124 ! Copy over "offline" settings 118 125 CALL init_phystokenc(offline,istphy) 126 127 ! Initialization for slab heat transport 128 type_ocean="force" 129 CALL getin_p('type_ocean',type_ocean) 130 slab_hdiff=.FALSE. 131 CALL getin_p('slab_hdiff',slab_hdiff) 132 slab_ekman=0 133 CALL getin_p('slab_ekman',slab_ekman) 134 IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman.GT.0))) THEN 135 CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,& 136 cu,cuvsurcv,cv,cvusurcu, & 137 aire,apoln,apols, & 138 aireu,airev) 139 END IF 119 140 120 141 ! Initialize tracer names, numbers, etc. for physics … … 159 180 #ifdef INCA 160 181 CALL init_inca_dim(klon_omp,nbp_lev,nbp_lon,nbp_lat - 1, & 161 rlonu ,rlatu,rlonv,rlatv)182 rlonudyn,rlatudyn,rlonvdyn,rlatvdyn) 162 183 #endif 163 184 END IF -
LMDZ5/branches/testing/libf/phydev/physiq_mod.F90
r2435 r2669 25 25 #ifdef CPP_XIOS 26 26 USE xios, ONLY: xios_update_calendar 27 USE wxios, only: wxios_add_vaxis, wxios_set_ timestep, wxios_closedef, &28 27 USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef 28 USE iophy, ONLY: histwrite_phy 29 29 #endif 30 30 … … 129 129 !XIOS 130 130 ! Declare available vertical axes to be used in output files: 131 !CALL wxios_add_vaxis("presnivs", "dummy-not-used", klev, presnivs)132 131 CALL wxios_add_vaxis("presnivs", klev, presnivs) 133 132 134 ! Declare time step length (in s):135 CALL wxios_set_ timestep(dtime)136 133 ! Declare calendar and time step 134 CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0) 135 137 136 !Finalize the context: 138 137 CALL wxios_closedef() 139 138 #endif 140 139 !$OMP END MASTER 140 !$OMP BARRIER 141 141 endif ! of if (debut) 142 142 -
LMDZ5/branches/testing/libf/phylmd/Dust/phytracr_spl_mod.F90
r2641 r2669 1170 1170 !print *,nbtr 1171 1171 do it=1,nbtr 1172 print *, it, tname(it+ 2)1173 if (tname(it+ 2) == 'PREC' ) then1172 print *, it, tname(it+nqo) 1173 if (tname(it+nqo) == 'PREC' ) then 1174 1174 id_prec=it 1175 1175 endif 1176 if (tname(it+ 2) == 'FINE' ) then1176 if (tname(it+nqo) == 'FINE' ) then 1177 1177 id_fine=it 1178 1178 endif 1179 if (tname(it+ 2) == 'COSS' ) then1179 if (tname(it+nqo) == 'COSS' ) then 1180 1180 id_coss=it 1181 1181 endif 1182 if (tname(it+ 2) == 'CODU' ) then1182 if (tname(it+nqo) == 'CODU' ) then 1183 1183 id_codu=it 1184 1184 endif 1185 if (tname(it+ 2) == 'SCDU' ) then1185 if (tname(it+nqo) == 'SCDU' ) then 1186 1186 id_scdu=it 1187 1187 endif … … 2621 2621 call iophys_ecrit('fav'//str2,1,'SOURCE','',source_tr(:,it)) 2622 2622 enddo 2623 #endif2624 2625 2623 do it=1,nbtr 2626 2624 write(str2,'(i2.2)') it 2627 2625 call iophys_ecrit('TRB'//str2,klev,'SOURCE','',tr_seri(:,:,it)) 2628 2626 enddo 2627 #endif 2629 2628 2630 2629 … … 2650 2649 2651 2650 #ifdef IOPHYS_DUST 2652 2653 2651 do it=1,nbtr 2654 2652 write(str2,'(i2.2)') it … … 2676 2674 ! 2677 2675 #ifdef IOPHYS_DUST 2678 !2679 2676 print *,'INPUT TO PRECUREMISSION' 2680 2677 call iophys_ecrit('ftsol',4,'ftsol','',ftsol) -
LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90
r2408 r2669 87 87 CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified) 88 88 CASE ('slab') 89 IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN 90 ! Read fraction from limit.nc 91 CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified) 92 ELSE 89 93 ! Get fraction from slab module 90 CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified) 94 CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified) 95 ENDIF 91 96 CASE ('couple') 92 97 ! Get fraction from the coupler -
LMDZ5/branches/testing/libf/phylmd/cloudth.F90
r2594 r2669 9 9 10 10 !=========================================================================== 11 ! Aut eur : Arnaud Octavio Jam (LMD/CNRS)11 ! Author : Arnaud Octavio Jam (LMD/CNRS) 12 12 ! Date : 25 Mai 2010 13 13 ! Objet : calcule les valeurs de qc et rneb dans les thermiques … … 49 49 REAL rneb(ngrid,klev) 50 50 REAL t(ngrid,klev) 51 REAL qsatmmussig1,qsatmmussig2,sqrt2pi, pi51 REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi 52 52 REAL rdd,cppd,Lv 53 53 REAL alth,alenv,ath,aenv 54 REAL sth,senv,sigma1s,sigma2s,xth,xenv 54 REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2 55 55 REAL Tbef,zdelta,qsatbef,zcor 56 56 REAL qlbef 57 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur 58 57 REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution 59 58 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) 60 59 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) … … 63 62 64 63 65 66 67 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!68 ! Gestion de deux versions de cloudth69 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!70 64 71 65 IF (iflag_cloudth_vert.GE.1) THEN … … 82 76 ! Initialisation des variables r?elles 83 77 !------------------------------------------------------------------------------- 84 sigma1(:,:)=0.85 sigma2(:,:)=0.86 qlth(:,:)=0.87 qlenv(:,:)=0.88 qltot(:,:)=0.89 rneb(:,:)=0.90 qcloud(:)=0.91 cth(:,:)=0.92 cenv(:,:)=0.93 ctot(:,:)=0.94 qsatmmussig1=0.95 qsatmmussig2=0.96 rdd=287.0497 cppd=1005.798 pi=3.1415999 Lv=2.5e6100 sqrt2pi=sqrt(2.*pi)101 102 103 104 !-------------------------------------------------------------------------------105 ! Calcul de la fraction du thermique et des ?cart-types des distributions106 !-------------------------------------------------------------------------------107 do ind1=1,ngrid108 109 if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then110 111 zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))112 113 114 ! zqenv(ind1)=po(ind1)115 Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)116 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))117 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)118 qsatbef=MIN(0.5,qsatbef)119 zcor=1./(1.-retv*qsatbef)120 qsatbef=qsatbef*zcor121 zqsatenv(ind1,ind2)=qsatbef122 123 124 125 126 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)127 aenv=1./(1.+(alenv*Lv/cppd))128 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))129 130 131 132 133 Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)134 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))135 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)136 qsatbef=MIN(0.5,qsatbef)137 zcor=1./(1.-retv*qsatbef)138 qsatbef=qsatbef*zcor139 zqsatth(ind1,ind2)=qsatbef140 141 alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)142 ath=1./(1.+(alth*Lv/cppd))143 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))144 145 146 147 !------------------------------------------------------------------------------148 ! Calcul des ?cart-types pour s149 !------------------------------------------------------------------------------150 151 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)152 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)153 ! if (paprs(ind1,ind2).gt.90000) then154 ! ratqs(ind1,ind2)=0.002155 ! else156 ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000157 ! endif158 sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)159 sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)160 ! sigma1s=ratqs(ind1,ind2)*po(ind1)161 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003162 163 !------------------------------------------------------------------------------164 ! Calcul de l'eau condens?e et de la couverture nuageuse165 !------------------------------------------------------------------------------166 sqrt2pi=sqrt(2.*pi)167 xth=sth/(sqrt(2.)*sigma2s)168 xenv=senv/(sqrt(2.)*sigma1s)169 cth(ind1,ind2)=0.5*(1.+1.*erf(xth))170 cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))171 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)172 173 qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))174 qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))175 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)176 177 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!178 if (ctot(ind1,ind2).lt.1.e-10) then179 ctot(ind1,ind2)=0.180 qcloud(ind1)=zqsatenv(ind1,ind2)181 182 else183 184 ctot(ind1,ind2)=ctot(ind1,ind2)185 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)186 187 endif188 189 190 ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'191 192 193 else ! gaussienne environnement seule194 195 zqenv(ind1)=po(ind1)196 Tbef=t(ind1,ind2)197 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))198 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)199 qsatbef=MIN(0.5,qsatbef)200 zcor=1./(1.-retv*qsatbef)201 qsatbef=qsatbef*zcor202 zqsatenv(ind1,ind2)=qsatbef203 204 205 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)206 zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)207 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)208 aenv=1./(1.+(alenv*Lv/cppd))209 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))210 211 212 sigma1s=ratqs(ind1,ind2)*zqenv(ind1)213 214 sqrt2pi=sqrt(2.*pi)215 xenv=senv/(sqrt(2.)*sigma1s)216 ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))217 qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))218 219 if (ctot(ind1,ind2).lt.1.e-3) then220 ctot(ind1,ind2)=0.221 qcloud(ind1)=zqsatenv(ind1,ind2)222 223 else224 225 ctot(ind1,ind2)=ctot(ind1,ind2)226 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)227 228 endif229 230 231 232 233 234 235 endif236 enddo237 238 return239 end240 241 242 243 !===========================================================================244 SUBROUTINE cloudth_vert(ngrid,klev,ind2, &245 & ztv,po,zqta,fraca, &246 & qcloud,ctot,zpspsk,paprs,ztla,zthl, &247 & ratqs,zqs,t)248 249 !===========================================================================250 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)251 ! Date : 25 Mai 2010252 ! Objet : calcule les valeurs de qc et rneb dans les thermiques253 !===========================================================================254 255 256 USE ioipsl_getin_p_mod, ONLY : getin_p257 258 IMPLICIT NONE259 260 #include "YOMCST.h"261 #include "YOETHF.h"262 #include "FCTTRE.h"263 #include "thermcell.h"264 #include "nuage.h"265 266 INTEGER itap,ind1,ind2267 INTEGER ngrid,klev,klon,l,ig268 269 REAL ztv(ngrid,klev)270 REAL po(ngrid)271 REAL zqenv(ngrid)272 REAL zqta(ngrid,klev)273 274 REAL fraca(ngrid,klev+1)275 REAL zpspsk(ngrid,klev)276 REAL paprs(ngrid,klev+1)277 REAL ztla(ngrid,klev)278 REAL zthl(ngrid,klev)279 280 REAL zqsatth(ngrid,klev)281 REAL zqsatenv(ngrid,klev)282 283 284 REAL sigma1(ngrid,klev)285 REAL sigma2(ngrid,klev)286 REAL qlth(ngrid,klev)287 REAL qlenv(ngrid,klev)288 REAL qltot(ngrid,klev)289 REAL cth(ngrid,klev)290 REAL cenv(ngrid,klev)291 REAL ctot(ngrid,klev)292 REAL rneb(ngrid,klev)293 REAL t(ngrid,klev)294 REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi295 REAL rdd,cppd,Lv,sqrt2,sqrtpi296 REAL alth,alenv,ath,aenv297 REAL sth,senv,sigma1s,sigma2s,xth,xenv298 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv299 REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth300 REAL Tbef,zdelta,qsatbef,zcor301 REAL qlbef302 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur303 ! Change the width of the PDF used for vertical subgrid scale heterogeneity304 ! (J Jouhaud, JL Dufresne, JB Madeleine)305 REAL,SAVE :: vert_alpha306 LOGICAL, SAVE :: firstcall = .TRUE.307 308 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)309 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)310 REAL zqs(ngrid), qcloud(ngrid)311 REAL erf312 313 !------------------------------------------------------------------------------314 ! Initialisation des variables r?elles315 !------------------------------------------------------------------------------316 78 sigma1(:,:)=0. 317 79 sigma2(:,:)=0. … … 334 96 sqrtpi=sqrt(pi) 335 97 336 IF (firstcall) THEN 337 vert_alpha=0.5 338 CALL getin_p('cloudth_vert_alpha',vert_alpha) 339 WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha 340 firstcall=.FALSE. 341 ENDIF 342 343 !------------------------------------------------------------------------------- 344 ! Calcul de la fraction du thermique et des ?cart-types des distributions 98 99 !------------------------------------------------------------------------------- 100 ! Cloud fraction in the thermals and standard deviation of the PDFs 345 101 !------------------------------------------------------------------------------- 346 102 do ind1=1,ngrid … … 350 106 zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) 351 107 352 353 ! zqenv(ind1)=po(ind1)354 108 Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) 355 109 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 356 qsatbef= R2ES *FOEEW(Tbef,zdelta)/paprs(ind1,ind2)110 qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) 357 111 qsatbef=MIN(0.5,qsatbef) 358 112 zcor=1./(1.-retv*qsatbef) … … 361 115 362 116 363 364 365 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 366 aenv=1./(1.+(alenv*Lv/cppd)) 367 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 368 369 117 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 118 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 119 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 120 121 !po = qt de l'environnement ET des thermique 122 !zqenv = qt environnement 123 !zqsatenv = qsat environnement 124 !zthl = Tl environnement 370 125 371 126 … … 378 133 zqsatth(ind1,ind2)=qsatbef 379 134 380 alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) 381 ath=1./(1.+(alth*Lv/cppd)) 382 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 383 384 385 386 !------------------------------------------------------------------------------ 387 ! Calcul des ?cart-types pour s 388 !------------------------------------------------------------------------------ 389 390 sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) 391 sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 135 alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 136 ath=1./(1.+(alth*Lv/cppd)) !al, p84 137 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 138 139 !zqta = qt thermals 140 !zqsatth = qsat thermals 141 !ztla = Tl thermals 142 143 !------------------------------------------------------------------------------ 144 ! s standard deviations 145 !------------------------------------------------------------------------------ 146 147 ! tests 148 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) 149 ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) 150 ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) 151 ! final option 152 sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) 153 sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 154 155 !------------------------------------------------------------------------------ 156 ! Condensed water and cloud cover 157 !------------------------------------------------------------------------------ 158 xth=sth/(sqrt2*sigma2s) 159 xenv=senv/(sqrt2*sigma1s) 160 cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam 161 cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam 162 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 163 164 qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) 165 qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) 166 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) 167 168 if (ctot(ind1,ind2).lt.1.e-10) then 169 ctot(ind1,ind2)=0. 170 qcloud(ind1)=zqsatenv(ind1,ind2) 171 else 172 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) 173 endif 174 175 else ! Environnement only, follow the if l.110 176 177 zqenv(ind1)=po(ind1) 178 Tbef=t(ind1,ind2) 179 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 180 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) 181 qsatbef=MIN(0.5,qsatbef) 182 zcor=1./(1.-retv*qsatbef) 183 qsatbef=qsatbef*zcor 184 zqsatenv(ind1,ind2)=qsatbef 185 186 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) 187 zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) 188 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 189 aenv=1./(1.+(alenv*Lv/cppd)) 190 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 191 192 sigma1s=ratqs(ind1,ind2)*zqenv(ind1) 193 194 xenv=senv/(sqrt2*sigma1s) 195 ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 196 qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) 197 198 if (ctot(ind1,ind2).lt.1.e-3) then 199 ctot(ind1,ind2)=0. 200 qcloud(ind1)=zqsatenv(ind1,ind2) 201 else 202 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) 203 endif 204 205 206 endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 207 enddo ! from the loop on ngrid l.108 208 return 209 end 210 211 212 213 !=========================================================================== 214 SUBROUTINE cloudth_vert(ngrid,klev,ind2, & 215 & ztv,po,zqta,fraca, & 216 & qcloud,ctot,zpspsk,paprs,ztla,zthl, & 217 & ratqs,zqs,t) 218 219 !=========================================================================== 220 ! Auteur : Arnaud Octavio Jam (LMD/CNRS) 221 ! Date : 25 Mai 2010 222 ! Objet : calcule les valeurs de qc et rneb dans les thermiques 223 !=========================================================================== 224 225 226 USE ioipsl_getin_p_mod, ONLY : getin_p 227 228 IMPLICIT NONE 229 230 #include "YOMCST.h" 231 #include "YOETHF.h" 232 #include "FCTTRE.h" 233 #include "thermcell.h" 234 #include "nuage.h" 235 236 INTEGER itap,ind1,ind2 237 INTEGER ngrid,klev,klon,l,ig 238 239 REAL ztv(ngrid,klev) 240 REAL po(ngrid) 241 REAL zqenv(ngrid) 242 REAL zqta(ngrid,klev) 243 244 REAL fraca(ngrid,klev+1) 245 REAL zpspsk(ngrid,klev) 246 REAL paprs(ngrid,klev+1) 247 REAL ztla(ngrid,klev) 248 REAL zthl(ngrid,klev) 249 250 REAL zqsatth(ngrid,klev) 251 REAL zqsatenv(ngrid,klev) 252 253 254 REAL sigma1(ngrid,klev) 255 REAL sigma2(ngrid,klev) 256 REAL qlth(ngrid,klev) 257 REAL qlenv(ngrid,klev) 258 REAL qltot(ngrid,klev) 259 REAL cth(ngrid,klev) 260 REAL cenv(ngrid,klev) 261 REAL ctot(ngrid,klev) 262 REAL rneb(ngrid,klev) 263 REAL t(ngrid,klev) 264 REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi 265 REAL rdd,cppd,Lv 266 REAL alth,alenv,ath,aenv 267 REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 268 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv 269 REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth 270 REAL Tbef,zdelta,qsatbef,zcor 271 REAL qlbef 272 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur 273 ! Change the width of the PDF used for vertical subgrid scale heterogeneity 274 ! (J Jouhaud, JL Dufresne, JB Madeleine) 275 REAL,SAVE :: vert_alpha 276 LOGICAL, SAVE :: firstcall = .TRUE. 277 278 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) 279 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 280 REAL zqs(ngrid), qcloud(ngrid) 281 REAL erf 282 283 !------------------------------------------------------------------------------ 284 ! Initialize 285 !------------------------------------------------------------------------------ 286 sigma1(:,:)=0. 287 sigma2(:,:)=0. 288 qlth(:,:)=0. 289 qlenv(:,:)=0. 290 qltot(:,:)=0. 291 rneb(:,:)=0. 292 qcloud(:)=0. 293 cth(:,:)=0. 294 cenv(:,:)=0. 295 ctot(:,:)=0. 296 qsatmmussig1=0. 297 qsatmmussig2=0. 298 rdd=287.04 299 cppd=1005.7 300 pi=3.14159 301 Lv=2.5e6 302 sqrt2pi=sqrt(2.*pi) 303 sqrt2=sqrt(2.) 304 sqrtpi=sqrt(pi) 305 306 IF (firstcall) THEN 307 vert_alpha=0.5 308 CALL getin_p('cloudth_vert_alpha',vert_alpha) 309 WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha 310 firstcall=.FALSE. 311 ENDIF 312 313 !------------------------------------------------------------------------------- 314 ! Calcul de la fraction du thermique et des ecart-types des distributions 315 !------------------------------------------------------------------------------- 316 do ind1=1,ngrid 317 318 if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement 319 320 zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv 321 322 323 Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) 324 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 325 qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) 326 qsatbef=MIN(0.5,qsatbef) 327 zcor=1./(1.-retv*qsatbef) 328 qsatbef=qsatbef*zcor 329 zqsatenv(ind1,ind2)=qsatbef 330 331 332 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 333 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 334 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 335 336 !zqenv = qt environnement 337 !zqsatenv = qsat environnement 338 !zthl = Tl environnement 339 340 341 Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) 342 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 343 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) 344 qsatbef=MIN(0.5,qsatbef) 345 zcor=1./(1.-retv*qsatbef) 346 qsatbef=qsatbef*zcor 347 zqsatth(ind1,ind2)=qsatbef 348 349 alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 350 ath=1./(1.+(alth*Lv/cppd)) !al, p84 351 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 352 353 354 !zqta = qt thermals 355 !zqsatth = qsat thermals 356 !ztla = Tl thermals 357 358 !------------------------------------------------------------------------------ 359 ! s standard deviation 360 !------------------------------------------------------------------------------ 361 362 sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) 363 sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) 364 ! tests 365 ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) 366 ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) 367 ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 368 ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) 392 369 ! if (paprs(ind1,ind2).gt.90000) then 393 370 ! ratqs(ind1,ind2)=0.002 … … 400 377 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 401 378 402 !------------------------------------------------------------------------------403 ! Calcul de l'eau condens?e et de la couverture nuageuse404 !------------------------------------------------------------------------------405 sqrt2pi=sqrt(2.*pi)406 xth=sth/(sqrt(2.)*sigma2s)407 xenv=senv/(sqrt(2.)*sigma1s)408 cth(ind1,ind2)=0.5*(1.+1.*erf(xth))409 cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))410 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)411 412 qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))413 qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))414 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)415 416 379 IF (iflag_cloudth_vert == 1) THEN 417 380 !------------------------------------------------------------------------------- 418 ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs 419 !------------------------------------------------------------------------------- 420 ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) 421 ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) 381 ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs 382 !------------------------------------------------------------------------------- 383 422 384 deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) 423 385 deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) 424 ! deltasenv=aenv*0.01*po(ind1) 425 ! deltasth=ath*0.01*zqta(ind1,ind2) 386 426 387 xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) 427 388 xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) … … 435 396 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 436 397 398 ! Environment 437 399 IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) 438 400 IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) … … 441 403 442 404 qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 443 ! qlenv(ind1,ind2)=IntJ 444 ! print*, qlenv(ind1,ind2),'VERIF EAU' 445 446 405 406 ! Thermal 447 407 IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) 448 ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))449 ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))450 408 IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) 451 409 IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) 452 410 IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) 453 411 qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 454 ! qlth(ind1,ind2)=IntJ455 ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'456 412 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) 457 413 … … 459 415 460 416 !------------------------------------------------------------------------------- 461 ! Version 3: Modification Jean Jouhaud. On condense a partir de-delta s417 ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s 462 418 !------------------------------------------------------------------------------- 463 419 ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) … … 470 426 xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) 471 427 xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) 428 exp_xenv1 = exp(-1.*xenv1**2) 429 exp_xenv2 = exp(-1.*xenv2**2) 472 430 xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) 473 431 xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) 474 ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)475 ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)432 exp_xth1 = exp(-1.*xth1**2) 433 exp_xth2 = exp(-1.*xth2**2) 476 434 477 435 cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) … … 479 437 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 480 438 481 IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2) 439 440 !environnement 441 IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 442 if (deltasenv .lt. 1.e-10) then 443 qlenv(ind1,ind2)=IntJ 444 else 482 445 IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) 483 IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) 484 IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2)) 485 486 ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) 487 ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) 488 ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) 489 446 IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) 447 IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) 490 448 qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 491 ! qlenv(ind1,ind2)=IntJ 492 ! print*, qlenv(ind1,ind2),'VERIF EAU' 493 494 IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2) 449 endif 450 451 452 !thermique 453 IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 454 if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!! 455 qlth(ind1,ind2)=IntJ 456 else 495 457 IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) 496 IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) 497 IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2)) 498 458 IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) 459 IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) 499 460 qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 500 ! qlth(ind1,ind2)=IntJ 501 ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' 461 endif 462 502 463 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) 503 504 505 464 506 465 507 466 ENDIF ! of if (iflag_cloudth_vert==1 or 2) 508 509 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!510 467 511 468 if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then … … 515 472 else 516 473 517 ctot(ind1,ind2)=ctot(ind1,ind2)518 474 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) 519 475 ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & … … 521 477 522 478 endif 523 524 525 526 ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' 527 528 529 else ! gaussienne environnement seule 479 480 else ! Environment only 530 481 531 482 zqenv(ind1)=po(ind1) … … 539 490 540 491 541 ! 492 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) 542 493 zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) 543 494 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) … … 548 499 sigma1s=ratqs(ind1,ind2)*zqenv(ind1) 549 500 550 sqrt2pi=sqrt(2.*pi) 551 xenv=senv/(sqrt(2.)*sigma1s) 501 xenv=senv/(sqrt2*sigma1s) 552 502 ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 553 qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt (2.)*cenv(ind1,ind2))503 qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) 554 504 555 505 if (ctot(ind1,ind2).lt.1.e-3) then … … 559 509 else 560 510 561 ctot(ind1,ind2)=ctot(ind1,ind2)562 511 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) 563 512 564 513 endif 565 514 566 567 568 569 570 571 endif 572 enddo 515 endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 516 enddo ! from the loop on ngrid l.333 573 517 574 518 return -
LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90
r2594 r2669 19 19 ok_ade, ok_aie, ok_cdnc, aerosol_couple, & 20 20 flag_aerosol, flag_aerosol_strat, new_aod, & 21 bl95_b0, bl95_b1,&21 flag_bc_internal_mixture, bl95_b0, bl95_b1,& 22 22 read_climoz, & 23 23 alp_offset) … … 64 64 ! ok_cdnc, ok cloud droplet number concentration 65 65 ! flag_aerosol_strat : flag pour les aerosols stratos 66 ! flag_bc_internal_mixture : use BC internal mixture if true 66 67 ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 67 68 ! … … 77 78 INTEGER :: flag_aerosol 78 79 INTEGER :: flag_aerosol_strat 80 LOGICAL :: flag_bc_internal_mixture 79 81 LOGICAL :: new_aod 80 82 REAL :: bl95_b0, bl95_b1 … … 95 97 INTEGER, SAVE :: flag_aerosol_omp 96 98 INTEGER, SAVE :: flag_aerosol_strat_omp 99 LOGICAL, SAVE :: flag_bc_internal_mixture_omp 97 100 LOGICAL, SAVE :: new_aod_omp 98 101 REAL,SAVE :: bl95_b0_omp, bl95_b1_omp … … 398 401 flag_aerosol_strat_omp = 0 399 402 CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp) 403 404 ! 405 !Config Key = flag_bc_internal_mixture 406 !Config Desc = state of mixture for BC aerosols 407 ! - n = external mixture 408 ! - y = internal mixture 409 !Config Def = n 410 !Config Help = Used in physiq.F / aeropt 411 ! 412 flag_bc_internal_mixture_omp = .false. 413 CALL getin('flag_bc_internal_mixture',flag_bc_internal_mixture_omp) 400 414 401 415 ! Temporary variable for testing purpose! … … 2136 2150 flag_aerosol=flag_aerosol_omp 2137 2151 flag_aerosol_strat=flag_aerosol_strat_omp 2152 flag_bc_internal_mixture=flag_bc_internal_mixture_omp 2138 2153 new_aod=new_aod_omp 2139 2154 aer_type = aer_type_omp … … 2311 2326 IF (ok_aie .AND. .NOT. ok_cdnc) THEN 2312 2327 CALL abort_physic('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1) 2328 ENDIF 2329 2330 ! BC internal mixture is only possible with RRTM & NSW=6 & flag_aerosol=6 or aerosol_couple 2331 IF (flag_bc_internal_mixture .AND. NSW.NE.6) THEN 2332 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with NSW=6',1) 2333 ENDIF 2334 IF (flag_bc_internal_mixture .AND. iflag_rrtm.NE.1) THEN 2335 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with RRTM',1) 2336 ENDIF 2337 IF (flag_bc_internal_mixture .AND. flag_aerosol.NE.6) THEN 2338 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with flag_aerosol=6',1) 2313 2339 ENDIF 2314 2340 -
LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90
r2594 r2669 42 42 "clmcalipso", "Lidar Mid-level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 43 43 TYPE(ctrl_outcosp), SAVE :: o_clhcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 44 "clhcalipso", "Lidar High t-level Cloud Fraction", "%", (/ ('', i=1, 3) /))44 "clhcalipso", "Lidar High-level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 45 45 TYPE(ctrl_outcosp), SAVE :: o_cltcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 46 46 "cltcalipso", "Lidar Total Cloud Fraction", "%", (/ ('', i=1, 3) /)) … … 145 145 "clmmodis", "MODIS Mid-level Cloud Fraction", "1", (/ ('', i=1, 3) /)) 146 146 TYPE(ctrl_outcosp), SAVE :: o_clhmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 147 "clhmodis", "MODIS High t-level Cloud Fraction", "1", (/ ('', i=1, 3) /))147 "clhmodis", "MODIS High-level Cloud Fraction", "1", (/ ('', i=1, 3) /)) 148 148 TYPE(ctrl_outcosp), SAVE :: o_cltmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 149 149 "cltmodis", "MODIS Total Cloud Fraction", "1", (/ ('', i=1, 3) /)) -
LMDZ5/branches/testing/libf/phylmd/cva_driver.F90
r2641 r2669 1077 1077 DO k = 1,nd 1078 1078 write (6, '(i4,5(1x,e13.6))'), & 1079 k, mp(i dcum(igout),k), water(idcum(igout),k), ice(idcum(igout),k), &1080 evap(i dcum(igout),k), fondue(idcum(igout),k)1079 k, mp(igout,k), water(igout,k), ice(igout,k), & 1080 evap(igout,k), fondue(igout,k) 1081 1081 ENDDO 1082 1082 Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainM ' 1083 1083 DO k = 1,nd 1084 1084 write (6, '(i4,2(1x,e13.6))'), & 1085 k, wdtrainA(i dcum(igout),k), wdtrainM(idcum(igout),k)1085 k, wdtrainA(igout,k), wdtrainM(igout,k) 1086 1086 ENDDO 1087 1087 ENDIF … … 1129 1129 ! 1130 1130 IF (debut) THEN 1131 PRINT *, ' cv3_yield -> fqd(1) = ', fqd(i dcum(igout), 1)1131 PRINT *, ' cv3_yield -> fqd(1) = ', fqd(igout, 1) 1132 1132 END IF !(debut) THEN 1133 1133 ! 1134 1134 IF (prt_level >= 10) THEN 1135 1135 Print *, 'cva_driver after cv3_yield:ft(1) , ftd(1) ', & 1136 ft(i dcum(igout),1), ftd(idcum(igout),1)1136 ft(igout,1), ftd(igout,1) 1137 1137 Print *, 'cva_driver after cv3_yield:fq(1) , fqd(1) ', & 1138 fq(i dcum(igout),1), fqd(idcum(igout),1)1138 fq(igout,1), fqd(igout,1) 1139 1139 ENDIF 1140 1140 ! -
LMDZ5/branches/testing/libf/phylmd/dimphy.F90
r2073 r2669 9 9 INTEGER,SAVE :: klevm1 10 10 INTEGER,SAVE :: kflev 11 INTEGER,SAVE :: nslay12 11 13 !$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon ,nslay)12 !$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon) 14 13 REAL,save,allocatable,dimension(:) :: zmasq 15 14 !$OMP THREADPRIVATE(zmasq) … … 24 23 25 24 klon=klon0 26 nslay=1 ! Slab, provisoire (F. Codron)27 25 kdlon=klon 28 26 kidia=1 -
LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90
r2542 r2669 624 624 zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i))) 625 625 else if (iflag_t_glace.ge.1) then 626 zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i))) 626 zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i))) 627 ! BUG corrige par JYG zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i))) 627 628 endif 628 629 endif -
LMDZ5/branches/testing/libf/phylmd/iophy.F90
r2542 r2669 17 17 #ifdef CPP_XIOS 18 18 INTERFACE histwrite_phy 19 MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios 19 !#ifdef CPP_XIOSnew 20 MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios,histwrite0d_xios 21 !#else 22 ! MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios 23 !#endif 24 20 25 END INTERFACE 21 26 #else … … 1331 1336 IF (prt_level >= 10) write(lunout,*)'End histrwrite3d_xios ',field_name 1332 1337 END SUBROUTINE histwrite3d_xios 1338 1339 #ifdef CPP_XIOS 1340 SUBROUTINE histwrite0d_xios(field_name, field) 1341 USE xios, only: xios_send_field 1342 IMPLICIT NONE 1343 1344 CHARACTER(LEN=*), INTENT(IN) :: field_name 1345 REAL, INTENT(IN) :: field ! --> scalar 1346 1347 !$OMP MASTER 1348 CALL xios_send_field(field_name, field) 1349 !$OMP END MASTER 1350 1351 END SUBROUTINE histwrite0d_xios 1352 #endif 1353 1333 1354 #endif 1334 1355 end module iophy -
LMDZ5/branches/testing/libf/phylmd/limit_slab.F90
r2542 r2669 164 164 diff_sst(:) = diff_sst_save(:) 165 165 diff_siv(:) = diff_siv_save(:) 166 ! For Debug purpose167 ! PRINT *,'limit_slab sst',MINVAL(diff_sst(:)),MAXVAL(diff_sst(:))168 ! PRINT *,'limit_slab siv',MINVAL(diff_siv(:)),MAXVAL(diff_siv(:))169 ! PRINT *,'limit_slab bils',MINVAL(lmt_bils(:)),MAXVAL(lmt_bils(:))170 166 171 167 END SUBROUTINE limit_slab -
LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90
r2408 r2669 9 9 USE indice_sol_mod 10 10 USE surface_data 11 USE mod_grid_phy_lmdz, ONLY: klon_glo 12 USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root 11 13 12 14 IMPLICIT NONE … … 14 16 PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice 15 17 16 !*********************************************************************************** *****18 !*********************************************************************************** 17 19 ! Global saved variables 18 !**************************************************************************************** 19 20 !*********************************************************************************** 21 ! number of slab vertical layers 22 INTEGER, PUBLIC, SAVE :: nslay 23 !$OMP THREADPRIVATE(nslay) 24 ! timestep for coupling (update slab temperature) in timesteps 20 25 INTEGER, PRIVATE, SAVE :: cpl_pas 21 26 !$OMP THREADPRIVATE(cpl_pas) 27 ! cyang = 1/heat capacity of top layer (rho.c.H) 22 28 REAL, PRIVATE, SAVE :: cyang 23 29 !$OMP THREADPRIVATE(cyang) 30 ! depth of slab layers (1 or 2) 24 31 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh 25 32 !$OMP THREADPRIVATE(slabh) 33 ! slab temperature 26 34 REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab 27 35 !$OMP THREADPRIVATE(tslab) 36 ! heat flux convergence due to Ekman 37 REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman 38 !$OMP THREADPRIVATE(dt_ekman) 39 ! heat flux convergence due to horiz diffusion 40 REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff 41 !$OMP THREADPRIVATE(dt_hdiff) 42 ! fraction of ocean covered by sea ice (sic / (oce+sic)) 28 43 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic 29 44 !$OMP THREADPRIVATE(fsic) 45 ! temperature of the sea ice 30 46 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice 31 47 !$OMP THREADPRIVATE(tice) 48 ! sea ice thickness, in kg/m2 32 49 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice 33 50 !$OMP THREADPRIVATE(seaice) 51 ! net surface heat flux, weighted by open ocean fraction 34 52 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bils 35 53 !$OMP THREADPRIVATE(slab_bils) 54 ! slab_bils accumulated over cpl_pas timesteps 36 55 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum 37 56 !$OMP THREADPRIVATE(bils_cum) 57 ! net heat flux into the ocean below the ice : conduction + solar radiation 38 58 REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg 39 59 !$OMP THREADPRIVATE(slab_bilg) 60 ! slab_bilg over cpl_pas timesteps 40 61 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum 41 62 !$OMP THREADPRIVATE(bilg_cum) 42 43 !**************************************************************************************** 63 ! wind stress saved over cpl_pas timesteps 64 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum 65 !$OMP THREADPRIVATE(taux_cum) 66 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum 67 !$OMP THREADPRIVATE(tauy_cum) 68 69 !*********************************************************************************** 44 70 ! Parameters (could be read in def file: move to slab_init) 45 !*********************************************************************************** *****71 !*********************************************************************************** 46 72 ! snow and ice physical characteristics: 47 73 REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp 48 74 REAL, PARAMETER :: t_melt=273.15 ! melting ice temp 49 75 REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3 50 REAL, PARAMETER :: ice_den=917. 51 REAL, PARAMETER :: sea_den=1025. 52 REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity 53 REAL, PARAMETER :: sno_cond=0.31*sno_den 76 REAL, PARAMETER :: ice_den=917. ! ice density 77 REAL, PARAMETER :: sea_den=1025. ! sea water density 78 REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice 79 REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow 54 80 REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice 81 REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, snow and ice 55 82 REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice 56 83 … … 74 101 REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice 75 102 REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice 76 REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow) 103 REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the 104 ! ice (no snow) 77 105 REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1) 78 106 79 !**************************************************************************************** 107 ! horizontal transport 108 LOGICAL, PUBLIC, SAVE :: slab_hdiff 109 !$OMP THREADPRIVATE(slab_hdiff) 110 REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion 111 !$OMP THREADPRIVATE(coef_hdiff) 112 INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment 113 !$OMP THREADPRIVATE(slab_ekman) 114 !$OMP THREADPRIVATE(slab_cadj) 115 116 !*********************************************************************************** 80 117 81 118 CONTAINS 82 119 ! 83 !*********************************************************************************** *****120 !*********************************************************************************** 84 121 ! 85 122 SUBROUTINE ocean_slab_init(dtime, pctsrf_rst) 86 123 !, seaice_rst etc 87 124 88 use IOIPSL 89 90 ! For ok_xxx vars (Ekman...) 91 INCLUDE "clesphys.h" 125 USE ioipsl_getin_p_mod, ONLY : getin_p 126 USE mod_phys_lmdz_transfert_para, ONLY : gather 127 USE slab_heat_transp_mod, ONLY : ini_slab_transp 92 128 93 129 ! Input variables 94 !*********************************************************************************** *****130 !*********************************************************************************** 95 131 REAL, INTENT(IN) :: dtime 96 132 ! Variables read from restart file 97 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst 133 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst 134 ! surface fractions from start file 98 135 99 136 ! Local variables 100 !************************************************************************************ ****137 !************************************************************************************ 101 138 INTEGER :: error 139 REAL, DIMENSION(klon_glo) :: zmasq_glo 102 140 CHARACTER (len = 80) :: abort_message 103 141 CHARACTER (len = 20) :: modname = 'ocean_slab_intit' 104 142 105 !**************************************************************************************** 143 !*********************************************************************************** 144 ! Define some parameters 145 !*********************************************************************************** 146 ! Number of slab layers 147 nslay=2 148 CALL getin_p('slab_layers',nslay) 149 print *,'number of slab layers : ',nslay 150 ! Layer thickness 151 ALLOCATE(slabh(nslay), stat = error) 152 IF (error /= 0) THEN 153 abort_message='Pb allocation slabh' 154 CALL abort_physic(modname,abort_message,1) 155 ENDIF 156 slabh(1)=50. 157 IF (nslay.GT.1) THEN 158 slabh(2)=150. 159 END IF 160 161 ! cyang = 1/heat capacity of top layer (rho.c.H) 162 cyang=1/(slabh(1)*sea_den*sea_cap) 163 164 ! cpl_pas coupling period (update of tslab and ice fraction) 165 ! pour un calcul a chaque pas de temps, cpl_pas=1 166 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour 167 CALL getin_p('cpl_pas',cpl_pas) 168 print *,'cpl_pas',cpl_pas 169 170 ! Horizontal diffusion 171 slab_hdiff=.FALSE. 172 CALL getin_p('slab_hdiff',slab_hdiff) 173 coef_hdiff=25000. 174 CALL getin_p('coef_hdiff',coef_hdiff) 175 ! Ekman transport 176 slab_ekman=0 177 CALL getin_p('slab_ekman',slab_ekman) 178 ! Convective adjustment 179 IF (nslay.EQ.1) THEN 180 slab_cadj=0 181 ELSE 182 slab_cadj=1 183 END IF 184 CALL getin_p('slab_cadj',slab_cadj) 185 186 !************************************************************************************ 106 187 ! Allocate surface fraction read from restart file 107 !************************************************************************************ ****188 !************************************************************************************ 108 189 ALLOCATE(fsic(klon), stat = error) 109 190 IF (error /= 0) THEN … … 112 193 ENDIF 113 194 fsic(:)=0. 195 !zmasq = continent fraction 114 196 WHERE (1.-zmasq(:)>EPSFRA) 115 197 fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:)) 116 198 END WHERE 117 199 118 !************************************************************************************ ****119 ! Allocate saved variables120 !************************************************************************************ ****200 !************************************************************************************ 201 ! Allocate saved fields 202 !************************************************************************************ 121 203 ALLOCATE(tslab(klon,nslay), stat=error) 122 204 IF (error /= 0) CALL abort_physic & … … 136 218 bils_cum(:) = 0.0 137 219 138 IF (version_ocean=='sicINT') THEN 220 IF (version_ocean=='sicINT') THEN ! interactive sea ice 139 221 ALLOCATE(slab_bilg(klon), stat = error) 140 222 IF (error /= 0) THEN … … 161 243 END IF 162 244 163 !**************************************************************************************** 164 ! Define some parameters 165 !**************************************************************************************** 166 ! Layer thickness 167 ALLOCATE(slabh(nslay), stat = error) 168 IF (error /= 0) THEN 169 abort_message='Pb allocation slabh' 170 CALL abort_physic(modname,abort_message,1) 245 IF (slab_hdiff) THEN !horizontal diffusion 246 ALLOCATE(dt_hdiff(klon,nslay), stat = error) 247 IF (error /= 0) THEN 248 abort_message='Pb allocation dt_hdiff' 249 CALL abort_physic(modname,abort_message,1) 250 ENDIF 251 dt_hdiff(:,:) = 0.0 171 252 ENDIF 172 slabh(1)=50. 173 ! cyang = 1/heat capacity of top layer (rho.c.H) 174 cyang=1/(slabh(1)*4.228e+06) 175 176 ! cpl_pas periode de couplage avec slab (update tslab, pctsrf) 177 ! pour un calcul à chaque pas de temps, cpl_pas=1 178 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour 179 CALL getin('cpl_pas',cpl_pas) 180 print *,'cpl_pas',cpl_pas 253 254 IF (slab_ekman.GT.0) THEN ! ekman transport 255 ALLOCATE(dt_ekman(klon,nslay), stat = error) 256 IF (error /= 0) THEN 257 abort_message='Pb allocation dt_ekman' 258 CALL abort_physic(modname,abort_message,1) 259 ENDIF 260 dt_ekman(:,:) = 0.0 261 ALLOCATE(taux_cum(klon), stat = error) 262 IF (error /= 0) THEN 263 abort_message='Pb allocation taux_cum' 264 CALL abort_physic(modname,abort_message,1) 265 ENDIF 266 taux_cum(:) = 0.0 267 ALLOCATE(tauy_cum(klon), stat = error) 268 IF (error /= 0) THEN 269 abort_message='Pb allocation tauy_cum' 270 CALL abort_physic(modname,abort_message,1) 271 ENDIF 272 tauy_cum(:) = 0.0 273 ENDIF 274 275 ! Initialize transport 276 IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN 277 CALL gather(zmasq,zmasq_glo) 278 !$OMP MASTER ! Only master thread 279 IF (is_mpi_root) THEN 280 CALL ini_slab_transp(zmasq_glo) 281 END IF 282 !$OMP END MASTER 283 END IF 181 284 182 285 END SUBROUTINE ocean_slab_init 183 286 ! 184 !*********************************************************************************** *****287 !*********************************************************************************** 185 288 ! 186 289 SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified) 187 290 188 USE limit_read_mod 189 190 ! INCLUDE "clesphys.h" 291 ! this routine sends back the sea ice and ocean fraction to the main physics 292 ! routine. Called only with interactive sea ice 191 293 192 294 ! Arguments 193 !************************************************************************************ ****194 INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant195 INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee196 REAL , INTENT(IN) :: dtime ! p as de temps de la physique (ens)295 !************************************************************************************ 296 INTEGER, INTENT(IN) :: itime ! current timestep 297 INTEGER, INTENT(IN) :: jour ! day in year (not 298 REAL , INTENT(IN) :: dtime ! physics timestep (s) 197 299 REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction 198 LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step 199 200 ! Local variables 201 !**************************************************************************************** 202 203 IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN 204 CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified) 205 ELSE 300 LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is 301 ! modified at this time step 302 206 303 pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:)) 207 304 pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:)) 208 305 is_modified=.TRUE. 209 END IF210 306 211 307 END SUBROUTINE ocean_slab_frac 212 308 ! 213 !************************************************************************************ ****309 !************************************************************************************ 214 310 ! 215 311 SUBROUTINE ocean_slab_noice( & … … 224 320 225 321 USE calcul_fluxs_mod 322 USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2 323 USE mod_phys_lmdz_para 226 324 227 325 INCLUDE "clesphys.h" 228 326 327 ! This routine 328 ! (1) computes surface turbulent fluxes over points with some open ocean 329 ! (2) reads additional Q-flux (everywhere) 330 ! (3) computes horizontal transport (diffusion & Ekman) 331 ! (4) updates slab temperature every cpl_pas ; creates new ice if needed. 332 333 ! Note : 334 ! klon total number of points 335 ! knon number of points with open ocean (varies with time) 336 ! knindex gives position of the knon points within klon. 337 ! In general, local saved variables have klon values 338 ! variables exchanged with PBL module have knon. 339 229 340 ! Input arguments 230 !**************************************************************************************** 231 INTEGER, INTENT(IN) :: itime 232 INTEGER, INTENT(IN) :: jour 233 INTEGER, INTENT(IN) :: knon 234 INTEGER, DIMENSION(klon), INTENT(IN) :: knindex 235 REAL, INTENT(IN) :: dtime 236 REAL, DIMENSION(klon), INTENT(IN) :: p1lay 237 REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm 238 REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 239 REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum 240 REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 241 REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 242 REAL, DIMENSION(klon), INTENT(IN) :: ps 243 REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness 244 REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in 245 REAL, DIMENSION(klon), INTENT(INOUT) :: radsol 341 !*********************************************************************************** 342 INTEGER, INTENT(IN) :: itime ! current timestep INTEGER, 343 INTEGER, INTENT(IN) :: jour ! day in year (for Q-Flux) 344 INTEGER, INTENT(IN) :: knon ! number of points 345 INTEGER, DIMENSION(klon), INTENT(IN) :: knindex 346 REAL, INTENT(IN) :: dtime ! timestep (s) 347 REAL, DIMENSION(klon), INTENT(IN) :: p1lay 348 REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm 349 ! drag coefficients 350 REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 351 REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum ! near surface T, q 352 REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 353 REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 354 ! exchange coefficients for boundary layer scheme 355 REAL, DIMENSION(klon), INTENT(IN) :: ps ! surface pressure 356 REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! surface wind 357 REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature 358 REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux 246 359 247 360 ! In/Output arguments 248 !************************************************************************************ ****249 REAL, DIMENSION(klon), INTENT(INOUT) :: snow 361 !************************************************************************************ 362 REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2 250 363 251 364 ! Output arguments 252 !************************************************************************************ ****365 !************************************************************************************ 253 366 REAL, DIMENSION(klon), INTENT(OUT) :: qsurf 254 367 REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat 255 368 REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 256 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new 369 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture 257 370 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l 258 371 REAL, DIMENSION(klon), INTENT(OUT) :: qflux 259 372 260 373 ! Local variables 261 !**************************************************************************************** 262 INTEGER :: i,ki 263 ! surface fluxes 374 !************************************************************************************ 375 INTEGER :: i,ki,k 376 REAL :: t_cadj 377 ! for surface heat fluxes 264 378 REAL, DIMENSION(klon) :: cal, beta, dif_grnd 265 ! f lux correction379 ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes 266 380 REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils 267 ! surface wind stress381 ! for surface wind stress 268 382 REAL, DIMENSION(klon) :: u0, v0 269 383 REAL, DIMENSION(klon) :: u1_lay, v1_lay 270 ! ice creation384 ! for new ice creation 271 385 REAL :: e_freeze, h_new, dfsic 272 273 !**************************************************************************************** 274 ! 1) Flux calculation 275 ! 276 !**************************************************************************************** 277 cal(:) = 0. 278 beta(:) = 1. 279 dif_grnd(:) = 0. 386 ! horizontal diffusion and Ekman local vars 387 ! dimension = global domain (klon_glo) instead of // subdomains 388 REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo 389 ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop 390 REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp 391 REAL, DIMENSION(klon_glo,nslay) :: tslab_glo 392 REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo 393 394 !**************************************************************************************** 395 ! 1) Surface fluxes calculation 396 ! 397 !**************************************************************************************** 398 cal(:) = 0. ! infinite thermal inertia 399 beta(:) = 1. ! wet surface 400 dif_grnd(:) = 0. ! no diffusion into ground 280 401 281 402 ! Suppose zero surface speed … … 285 406 v1_lay(:) = v1(:) - v0(:) 286 407 408 ! Compute latent & sensible fluxes 287 409 CALL calcul_fluxs(knon, is_oce, dtime, & 288 410 tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, & … … 292 414 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) 293 415 294 ! - Flux calculation at first modele level for U and V 295 CALL calcul_flux_wind(knon, dtime, & 296 u0, v0, u1, v1, gustiness, cdragm, & 297 AcoefU, AcoefV, BcoefU, BcoefV, & 298 p1lay, temp_air, & 299 flux_u1, flux_v1) 300 301 ! Accumulate total fluxes locally 416 ! save total cumulated heat fluxes locally 417 ! radiative + turbulent + melt of falling snow 302 418 slab_bils(:)=0. 303 419 DO i=1,knon … … 306 422 -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki))) 307 423 bils_cum(ki)=bils_cum(ki)+slab_bils(ki) 308 ! Also taux, tauy, saved vars...309 424 END DO 310 425 311 !**************************************************************************************** 312 ! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc 426 ! Compute surface wind stress 427 CALL calcul_flux_wind(knon, dtime, & 428 u0, v0, u1, v1, gustiness, cdragm, & 429 AcoefU, AcoefV, BcoefU, BcoefV, & 430 p1lay, temp_air, & 431 flux_u1, flux_v1) 432 433 ! save cumulated wind stress 434 IF (slab_ekman.GT.0) THEN 435 DO i=1,knon 436 ki=knindex(i) 437 taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas 438 tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas 439 END DO 440 ENDIF 441 442 !**************************************************************************************** 443 ! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc 313 444 ! 314 445 !**************************************************************************************** 315 446 lmt_bils(:)=0. 316 CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus447 CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) 317 448 ! lmt_bils and diff_sst,siv saved by limit_slab 318 449 qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400. … … 320 451 321 452 !**************************************************************************************** 322 ! 3) Recalculate new temperature 323 ! 453 ! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport) 454 ! Bring to freezing temp and make sea ice if necessary 455 ! 324 456 !***********************************************o***************************************** 325 457 tsurf_new=tsurf_in 326 458 IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction 327 ! Compute transport 328 ! Add read QFlux and SST tendency 329 tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 330 ! Add cumulated surface fluxes 331 tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime 332 ! Update surface temperature 333 SELECT CASE(version_ocean) 334 CASE('sicNO') 335 DO i=1,knon 336 ki=knindex(i) 337 tsurf_new(i)=tslab(ki,1) 459 ! *********************************** 460 ! Horizontal transport 461 ! *********************************** 462 IF (slab_ekman.GT.0) THEN 463 ! copy wind stress to global var 464 CALL gather(taux_cum,taux_glo) 465 CALL gather(tauy_cum,tauy_glo) 466 END IF 467 468 IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN 469 CALL gather(tslab,tslab_glo) 470 ! Compute horiz transport on one process only 471 !$OMP MASTER ! Only master thread 472 IF (is_mpi_root) THEN ! Only master processus 473 IF (slab_hdiff) THEN 474 dt_hdiff_glo(:,:)=0. 475 END IF 476 IF (slab_ekman.GT.0) THEN 477 dt_ekman_glo(:,:)=0. 478 END IF 479 DO i=1,cpl_pas ! time splitting for numerical stability 480 IF (slab_ekman.GT.0) THEN 481 SELECT CASE (slab_ekman) 482 CASE (1) 483 CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) 484 CASE (2) 485 CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) 486 CASE DEFAULT 487 dt_ekman_tmp(:,:)=0. 488 END SELECT 489 dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:) 490 ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1 491 DO k=1,nslay 492 dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den) 493 ENDDO 494 tslab_glo=tslab_glo+dt_ekman_tmp*dtime 495 ENDIF 496 IF (slab_hdiff) THEN ! horizontal diffusion 497 ! laplacian of slab T 498 CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp) 499 ! multiply by diff coef and normalize to 50m slab equivalent 500 dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh) 501 dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:) 502 tslab_glo=tslab_glo+dt_hdiff_tmp*dtime 503 END IF 504 END DO ! time splitting 505 IF (slab_hdiff) THEN 506 !dt_hdiff_glo saved in W/m2 507 DO k=1,nslay 508 dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas 338 509 END DO 339 CASE('sicOBS') ! check for sea ice or tslab below freezing 340 DO i=1,knon 341 ki=knindex(i) 342 IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN 343 tslab(ki,1)=t_freeze 344 END IF 345 tsurf_new(i)=tslab(ki,1) 346 END DO 347 CASE('sicINT') 348 DO i=1,knon 349 ki=knindex(i) 350 IF (fsic(ki).LT.epsfra) THEN ! Free of ice 351 IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice 352 ! quantity of new ice formed 353 e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat 354 ! new ice 355 tice(ki)=t_freeze 356 fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin) 357 IF (fsic(ki).GT.ice_frac_min) THEN 358 seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max) 359 tslab(ki,1)=t_freeze 360 ELSE 361 fsic(ki)=0. 362 END IF 363 tsurf_new(i)=t_freeze 364 ELSE 365 tsurf_new(i)=tslab(ki,1) 366 END IF 367 ELSE ! ice present 368 tsurf_new(i)=t_freeze 369 IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice 370 ! quantity of new ice formed over open ocean 371 e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) & 372 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) 373 ! new ice height and fraction 374 h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new 375 dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new) 376 h_new=MIN(e_freeze/dfsic,h_ice_max) 377 ! update tslab to freezing over open ocean only 378 tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki)) 379 ! update sea ice 380 seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) & 381 /(dfsic+fsic(ki)) 382 fsic(ki)=fsic(ki)+dfsic 383 ! update snow? 384 END IF !freezing 385 END IF ! sea ice present 386 END DO 387 END SELECT 388 bils_cum(:)=0.0! clear cumulated fluxes 510 END IF 511 IF (slab_ekman.GT.0) THEN 512 ! dt_ekman_glo saved in W/m2 513 dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas 514 END IF 515 END IF ! is_mpi_root 516 !$OMP END MASTER 517 !$OMP BARRIER 518 ! Send new fields back to all processes 519 CALL Scatter(tslab_glo,tslab) 520 IF (slab_hdiff) THEN 521 CALL Scatter(dt_hdiff_glo,dt_hdiff) 522 END IF 523 IF (slab_ekman.GT.0) THEN 524 CALL Scatter(dt_ekman_glo,dt_ekman) 525 ! clear wind stress 526 taux_cum(:)=0. 527 tauy_cum(:)=0. 528 END IF 529 ENDIF ! transport 530 531 ! *********************************** 532 ! Other heat fluxes 533 ! *********************************** 534 ! Add read QFlux 535 tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 536 ! Add cumulated surface fluxes 537 tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime 538 ! Convective adjustment if 2 layers 539 IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN 540 DO i=1,klon 541 IF (tslab(i,2).GT.tslab(i,1)) THEN 542 ! mean (mass-weighted) temperature 543 t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:)) 544 tslab(i,1)=t_cadj 545 tslab(i,2)=t_cadj 546 END IF 547 END DO 548 END IF 549 ! *********************************** 550 ! Update surface temperature and ice 551 ! *********************************** 552 SELECT CASE(version_ocean) 553 CASE('sicNO') ! no sea ice even below freezing ! 554 DO i=1,knon 555 ki=knindex(i) 556 tsurf_new(i)=tslab(ki,1) 557 END DO 558 CASE('sicOBS') ! "realistic" case, for prescribed sea ice 559 ! tslab cannot be below freezing, or above it if there is sea ice 560 DO i=1,knon 561 ki=knindex(i) 562 IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN 563 tslab(ki,1)=t_freeze 564 END IF 565 tsurf_new(i)=tslab(ki,1) 566 END DO 567 CASE('sicINT') ! interactive sea ice 568 DO i=1,knon 569 ki=knindex(i) 570 IF (fsic(ki).LT.epsfra) THEN ! Free of ice 571 IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice 572 ! quantity of new ice formed 573 e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat 574 ! new ice 575 tice(ki)=t_freeze 576 fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin) 577 IF (fsic(ki).GT.ice_frac_min) THEN 578 seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max) 579 tslab(ki,1)=t_freeze 580 ELSE 581 fsic(ki)=0. 582 END IF 583 tsurf_new(i)=t_freeze 584 ELSE 585 tsurf_new(i)=tslab(ki,1) 586 END IF 587 ELSE ! ice present 588 tsurf_new(i)=t_freeze 589 IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice 590 ! quantity of new ice formed over open ocean 591 e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) & 592 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) 593 ! new ice height and fraction 594 h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new 595 dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new) 596 h_new=MIN(e_freeze/dfsic,h_ice_max) 597 ! update tslab to freezing over open ocean only 598 tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki)) 599 ! update sea ice 600 seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) & 601 /(dfsic+fsic(ki)) 602 fsic(ki)=fsic(ki)+dfsic 603 ! update snow? 604 END IF ! tslab below freezing 605 END IF ! sea ice present 606 END DO 607 END SELECT 608 bils_cum(:)=0.0! clear cumulated fluxes 389 609 END IF ! coupling time 390 610 END SUBROUTINE ocean_slab_noice … … 530 750 seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime) 531 751 ENDIF 532 ! Melt / Freeze from above if Tsurf>0752 ! Melt / Freeze snow from above if Tsurf>0 533 753 IF (tsurf_new(i).GT.t_melt) THEN 534 ! energy available for melting snow (in kg /m2 of snow)754 ! energy available for melting snow (in kg of melted snow /m2) 535 755 e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & 536 756 /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) … … 672 892 END IF ! coupling time 673 893 674 !tests fraction glace894 !tests ice fraction 675 895 WHERE (fsic.LT.ice_frac_min) 676 896 tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang … … 691 911 IF (ALLOCATED(tslab)) DEALLOCATE(tslab) 692 912 IF (ALLOCATED(fsic)) DEALLOCATE(fsic) 913 IF (ALLOCATED(tice)) DEALLOCATE(tice) 914 IF (ALLOCATED(seaice)) DEALLOCATE(seaice) 693 915 IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils) 694 916 IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) 695 917 IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) 696 918 IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) 697 IF (ALLOCATED(tslab)) DEALLOCATE(tslab) 919 IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) 920 IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) 921 IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) 922 IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) 698 923 699 924 END SUBROUTINE ocean_slab_final -
LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90
r2471 r2669 3095 3095 endif 3096 3096 mfois(nsrf) = mfois(nsrf) + 1 3097 ! F. Codron sensible default values for ocean and sea ice 3098 IF (nsrf.EQ.is_oce) THEN 3099 tsurf(i,nsrf) = 271.35 ! Freezing sea water 3100 DO k=1,nsw 3101 alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo 3102 alb_dif(i,k,nsrf) = 0.06 3103 END DO 3104 ELSE IF (nsrf.EQ.is_sic) THEN 3105 tsurf(i,nsrf) = 273.15 ! Melting ice 3106 DO k=1,nsw 3107 alb_dir(i,k,nsrf) = 0.3 ! thin ice 3108 alb_dif(i,k,nsrf) = 0.3 3109 END DO 3110 END IF 3111 ! F. Codron 3097 3112 ELSE 3098 3113 ! The continents have changed. The new fraction receives the mean sum of the existent fractions -
LMDZ5/branches/testing/libf/phylmd/phyetat0.F90
r2641 r2669 3 3 SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0) 4 4 5 USE dimphy, only: klon, zmasq, klev , nslay5 USE dimphy, only: klon, zmasq, klev 6 6 USE iophy, ONLY : init_iophy_new 7 7 USE ocean_cpl_mod, ONLY : ocean_cpl_init … … 25 25 USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send 26 26 USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic 27 USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init27 USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init 28 28 USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy 29 29 … … 443 443 IF ( type_ocean == 'slab' ) THEN 444 444 CALL ocean_slab_init(dtime, pctsrf) 445 found=phyetat0_get(nslay,tslab,"tslab","tslab",0.) 445 IF (nslay.EQ.1) THEN 446 found=phyetat0_get(1,tslab,"tslab01","tslab",0.) 447 IF (.NOT. found) THEN 448 found=phyetat0_get(1,tslab,"tslab","tslab",0.) 449 END IF 450 ELSE 451 DO i=1,nslay 452 WRITE(str2,'(i2.2)') i 453 found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.) 454 END DO 455 END IF 446 456 IF (.NOT. found) THEN 447 457 PRINT*, "phyetat0: Le champ <tslab> est absent" … … 453 463 454 464 ! Sea ice variables 455 found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)456 465 IF (version_ocean == 'sicINT') THEN 466 found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.) 457 467 IF (.NOT. found) THEN 458 468 PRINT*, "phyetat0: Le champ <tice> est absent" … … 460 470 tice(:)=ftsol(:,is_sic) 461 471 END IF 472 found=phyetat0_get(1,seaice,"seaice","seaice",0.) 462 473 IF (.NOT. found) THEN 463 474 PRINT*, "phyetat0: Le champ <seaice> est absent" -
LMDZ5/branches/testing/libf/phylmd/phyredem.F90
r2641 r2669 30 30 USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra 31 31 USE surface_data, ONLY: type_ocean, version_ocean 32 USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic32 USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic 33 33 USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys 34 34 … … 58 58 59 59 INTEGER isoil, nsrf,isw 60 CHARACTER (len= 7) :: str760 CHARACTER (len=2) :: str2 61 61 CHARACTER (len=256) :: nam, lnam 62 62 INTEGER :: it, iiq … … 304 304 ! Restart variables for Slab ocean 305 305 IF (type_ocean == 'slab') THEN 306 CALL put_field("tslab", "Slab ocean temperature", tslab) 306 IF (nslay.EQ.1) THEN 307 CALL put_field("tslab", "Slab ocean temperature", tslab) 308 ELSE 309 DO it=1,nslay 310 WRITE(str2,'(i2.2)') it 311 CALL put_field("tslab"//str2, "Slab ocean temperature", tslab(:,it)) 312 END DO 313 END IF 307 314 IF (version_ocean == 'sicINT') THEN 308 315 CALL put_field("seaice", "Slab seaice (kg/m2)", seaice) -
LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90
r2641 r2669 677 677 TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), & 678 678 'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /)) 679 TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), & 680 'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 9) /)) 681 TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), & 682 'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 9) /)) 679 683 TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), & 680 684 'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /)) -
LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90
r2641 r2669 24 24 ! defined and initialised in phys_output_mod.F90 25 25 26 USE dimphy, only: klon, klev, klevp1 , nslay26 USE dimphy, only: klon, klev, klevp1 27 27 USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy 28 28 USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat … … 89 89 o_slab_qflux, o_tslab, o_slab_bils, & 90 90 o_slab_bilg, o_slab_sic, o_slab_tice, & 91 o_slab_hdiff, o_slab_ekman, & 91 92 o_weakinv, o_dthmin, o_cldtau, & 92 93 o_cldemi, o_pr_con_l, o_pr_con_i, & … … 291 292 292 293 293 USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice 294 USE ocean_slab_mod, only: nslay, tslab, slab_bils, slab_bilg, tice, & 295 seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff 294 296 USE pbl_surface_mod, only: snow 295 297 USE indice_sol_mod, only: nbsrf … … 312 314 USE phys_cal_mod, only : mth_len 313 315 316 #ifdef CPP_RRTM 317 USE YOESW, ONLY : RSUN 318 #endif 314 319 315 320 IMPLICIT NONE … … 355 360 #endif 356 361 REAL, PARAMETER :: un_jour=86400. 362 INTEGER ISW 363 CHARACTER*1 ch1 357 364 358 365 ! On calcul le nouveau tau: … … 403 410 CALL histwrite_phy(o_contfracOR, pctsrf(:,is_ter)) 404 411 CALL histwrite_phy(o_aireTER, paire_ter) 412 ! 413 #ifdef CPP_XIOS 414 CALL histwrite_phy("R_ecc",R_ecc) 415 CALL histwrite_phy("R_peri",R_peri) 416 CALL histwrite_phy("R_incl",R_incl) 417 CALL histwrite_phy("solaire",solaire) 418 ! 419 #ifdef CPP_RRTM 420 IF (iflag_rrtm.EQ.1) THEN 421 DO ISW=1, NSW 422 WRITE(ch1,'(i1)') ISW 423 ! zx_tmp_0d=RSUN(ISW) 424 ! CALL histwrite_phy("rsun"//ch1,zx_tmp_0d) 425 CALL histwrite_phy("rsun"//ch1,RSUN(ISW)) 426 ENDDO 427 ENDIF 428 #endif 429 ! 430 CALL histwrite_phy("co2_ppm",co2_ppm) 431 CALL histwrite_phy("CH4_ppb",CH4_ppb) 432 CALL histwrite_phy("N2O_ppb",N2O_ppb) 433 CALL histwrite_phy("CFC11_ppt",CFC11_ppt) 434 CALL histwrite_phy("CFC12_ppt",CFC12_ppt) 435 ! 436 #endif 437 405 438 !!! Champs 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 406 439 ! Simulateur AIRS … … 951 984 CALL histwrite_phy(o_tslab, zx_tmp_fi2d) 952 985 ELSE 953 CALL histwrite_phy(o_tslab, tslab )986 CALL histwrite_phy(o_tslab, tslab(:,1:nslay)) 954 987 END IF 955 988 IF (version_ocean=='sicINT') THEN … … 957 990 CALL histwrite_phy(o_slab_tice, tice) 958 991 CALL histwrite_phy(o_slab_sic, seaice) 992 END IF 993 IF (slab_hdiff) THEN 994 IF (nslay.EQ.1) THEN 995 zx_tmp_fi2d(:)=dt_hdiff(:,1) 996 CALL histwrite_phy(o_slab_hdiff, zx_tmp_fi2d) 997 ELSE 998 CALL histwrite_phy(o_slab_hdiff, dt_hdiff(:,1:nslay)) 999 END IF 1000 END IF 1001 IF (slab_ekman.GT.0) THEN 1002 IF (nslay.EQ.1) THEN 1003 zx_tmp_fi2d(:)=dt_ekman(:,1) 1004 CALL histwrite_phy(o_slab_ekman, zx_tmp_fi2d) 1005 ELSE 1006 CALL histwrite_phy(o_slab_ekman, dt_ekman(:,1:nslay)) 1007 END IF 959 1008 END IF 960 1009 ENDIF !type_ocean == force/slab -
LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90
r2641 r2669 233 233 USE ioipsl_getin_p_mod, ONLY : getin_p 234 234 235 #ifndef CPP_XIOS 235 236 USE paramLMDZ_phy_mod 237 #endif 236 238 237 239 USE cmp_seri_mod … … 582 584 LOGICAL,SAVE :: ok_adjwk=.FALSE. 583 585 !$OMP THREADPRIVATE(ok_adjwk) 584 REAL,SAVE :: oliqmax=999. 585 !$OMP THREADPRIVATE(oliqmax )586 REAL,SAVE :: oliqmax=999.,oicemax=999. 587 !$OMP THREADPRIVATE(oliqmax,oicemax) 586 588 REAL, SAVE :: alp_offset 587 589 !$OMP THREADPRIVATE(alp_offset) … … 983 985 INTEGER, SAVE :: flag_aerosol 984 986 !$OMP THREADPRIVATE(flag_aerosol) 987 LOGICAL, SAVE :: flag_bc_internal_mixture 988 !$OMP THREADPRIVATE(flag_bc_internal_mixture) 985 989 LOGICAL, SAVE :: new_aod 986 990 !$OMP THREADPRIVATE(new_aod) … … 1011 1015 ! edges of pressure intervals for ozone climatologies, in Pa, in strictly 1012 1016 ! ascending order 1013 1014 integer, save:: co3i = 01015 ! time index in NetCDF file of current ozone fields1016 !$OMP THREADPRIVATE(co3i)1017 1017 1018 1018 integer ro3i … … 1140 1140 ok_ade, ok_aie, ok_cdnc, aerosol_couple, & 1141 1141 flag_aerosol, flag_aerosol_strat, new_aod, & 1142 bl95_b0, bl95_b1, &1142 flag_bc_internal_mixture, bl95_b0, bl95_b1, & 1143 1143 ! nv flags pour la convection et les 1144 1144 ! poches froides … … 1197 1197 CALL getin_p('ok_adjwk',ok_adjwk) 1198 1198 CALL getin_p('oliqmax',oliqmax) 1199 CALL getin_p('oicemax',oicemax) 1199 1200 CALL getin_p('ratqsp0',ratqsp0) 1200 1201 CALL getin_p('ratqsdp',ratqsdp) … … 1551 1552 WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3) 1552 1553 1554 #ifndef CPP_XIOS 1553 1555 CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM) 1556 #endif 1554 1557 1555 1558 #endif … … 1906 1909 ! Prescrire l'ozone et calculer l'albedo sur l'ocean. 1907 1910 ! 1908 if (read_climoz >= 1) then 1909 ! Ozone from a file 1910 ! Update required ozone index: 1911 ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1 1912 if (ro3i == 361) ro3i = 360 1913 ! (This should never occur, except perhaps because of roundup 1914 ! error. See documentation.) 1915 if (ro3i /= co3i) then 1916 ! Update ozone field: 1917 if (read_climoz == 1) then 1918 call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, & 1919 press_in_edg=press_climoz, paprs=paprs, v3=wo) 1920 else 1921 ! read_climoz == 2 1922 call regr_pr_av(ncid_climoz, (/"tro3 ", & 1923 "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, & 1924 paprs=paprs, v3=wo) 1925 end if 1926 ! Convert from mole fraction of ozone to column density of ozone in a 1927 ! cell, in kDU: 1928 forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd & 1929 * zmasse / dobson_u / 1e3 1930 ! (By regridding ozone values for LMDZ only once every 360th of 1931 ! year, we have already neglected the variation of pressure in one 1932 ! 360th of year. So do not recompute "wo" at each time step even if 1933 ! "zmasse" changes a little.) 1934 co3i = ro3i 1935 end if 1936 ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN 1911 ! Update ozone if day change 1912 IF (MOD(itap-1,lmt_pas) == 0) THEN 1913 IF (read_climoz == 0) THEN 1937 1914 ! Once per day, update ozone from Royer: 1938 1939 IF (solarlong0<-999.) then 1940 ! Generic case with evolvoing season 1941 zzz=real(days_elapsed+1) 1942 ELSE IF (abs(solarlong0-1000.)<1.e-4) then 1943 ! Particular case with annual mean insolation 1944 zzz=real(90) ! could be revisited 1945 IF (read_climoz/=-1) THEN 1946 abort_message ='read_climoz=-1 is recommended when ' & 1947 // 'solarlong0=1000.' 1948 CALL abort_physic (modname,abort_message,1) 1949 ENDIF 1950 ELSE 1915 IF (solarlong0<-999.) then 1916 ! Generic case with evolvoing season 1917 zzz=real(days_elapsed+1) 1918 ELSE IF (abs(solarlong0-1000.)<1.e-4) then 1919 ! Particular case with annual mean insolation 1920 zzz=real(90) ! could be revisited 1921 IF (read_climoz/=-1) THEN 1922 abort_message ='read_climoz=-1 is recommended when ' & 1923 // 'solarlong0=1000.' 1924 CALL abort_physic (modname,abort_message,1) 1925 ENDIF 1926 ELSE 1951 1927 ! Case where the season is imposed with solarlong0 1952 1928 zzz=real(90) ! could be revisited 1953 ENDIF 1954 wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz) 1929 ENDIF 1930 1931 wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz) 1932 ELSE 1933 ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1 1934 if (ro3i == 361) ro3i = 360 1935 if (read_climoz == 1) then 1936 call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, & 1937 press_in_edg=press_climoz, paprs=paprs, v3=wo) 1938 else 1939 ! read_climoz == 2 1940 call regr_pr_av(ncid_climoz, (/"tro3 ", & 1941 "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, & 1942 paprs=paprs, v3=wo) 1943 end if 1944 ! Convert from mole fraction of ozone to column density of ozone in a 1945 ! cell, in kDU: 1946 forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd & 1947 * zmasse / dobson_u / 1e3 1948 ! (By regridding ozone values for LMDZ only once every 360th of 1949 ! year, we have already neglected the variation of pressure in one 1950 ! 360th of year. So do not recompute "wo" at each time step even if 1951 ! "zmasse" changes a little.) 1952 1953 ENDIF 1955 1954 ENDIF 1956 1955 ! … … 3034 3033 'lsc',abortphy,flag_inhib_tend) 3035 3034 rain_num(:)=0. 3036 DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse3035 DO k = 1, klev 3037 3036 DO i = 1, klon 3038 3037 IF (ql_seri(i,k)>oliqmax) THEN … … 3042 3041 ENDDO 3043 3042 ENDDO 3043 IF (nqo==3) THEN 3044 DO k = 1, klev 3045 DO i = 1, klon 3046 IF (qs_seri(i,k)>oicemax) THEN 3047 rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys 3048 qs_seri(i,k)=oicemax 3049 ENDIF 3050 ENDDO 3051 ENDDO 3052 ENDIF 3044 3053 3045 3054 !--------------------------------------------------------------------------- … … 3388 3397 ! 3389 3398 CALL readaerosol_optic_rrtm( debut, aerosol_couple, & 3390 new_aod, flag_aerosol, itap, jD_cur-jD_ref, &3399 new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & 3391 3400 pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & 3392 3401 tr_seri, mass_solu_aero, mass_solu_aero_pi, & … … 4517 4526 4518 4527 4528 #ifndef CPP_XIOS 4519 4529 CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync) 4530 #endif 4520 4531 4521 4532 #endif -
LMDZ5/branches/testing/libf/phylmd/readaerosol_optic.F90
r2408 r2669 166 166 m_allaer_pi(:,:,id_CSSO4M_phy) = 0. ! CSSO4M pre-ind 167 167 m_allaer_pi(:,:,id_SSSSM_phy) = sssupco_pi(:,:) ! SSSSM pre-ind 168 m_allaer_pi(:,:,id_ ASSSM_phy) = sscoarse_pi(:,:) ! CSSSM pre-ind169 m_allaer_pi(:,:,id_ CIDUSTM_phy)= ssacu_pi(:,:) ! ASSSM pre-ind170 m_allaer_pi(:,:,id_ AIBCM_phy)= cidust_pi(:,:) ! CIDUSTM pre-ind168 m_allaer_pi(:,:,id_CSSSM_phy) = sscoarse_pi(:,:) ! CSSSM pre-ind 169 m_allaer_pi(:,:,id_ASSSM_phy) = ssacu_pi(:,:) ! ASSSM pre-ind 170 m_allaer_pi(:,:,id_CIDUSTM_phy) = cidust_pi(:,:) ! CIDUSTM pre-ind 171 171 m_allaer_pi(:,:,id_AIBCM_phy) = bcins_pi(:,:) ! AIBCM pre-ind 172 172 m_allaer_pi(:,:,id_AIPOMM_phy) = pomins_pi(:,:) ! AIPOMM pre-ind -
LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90
r2641 r2669 3 3 ! 4 4 5 SUBROUTINE AEROPT_5WV_RRTM(& 6 pdel, m_allaer, delt, & 7 RHcl, ai, flag_aerosol, & 8 pplay, t_seri, & 5 SUBROUTINE AEROPT_5WV_RRTM( & 6 pdel, m_allaer, & 7 RHcl, ai, flag_aerosol, & 8 flag_bc_internal_mixture, & 9 pplay, t_seri, & 9 10 tausum, tau ) 10 11 … … 54 55 ! 55 56 REAL, DIMENSION(klon,klev), INTENT(in) :: pdel 56 REAL, INTENT(in) :: delt57 57 REAL, DIMENSION(klon,klev,naero_tot), INTENT(in) :: m_allaer 58 58 REAL, DIMENSION(klon,klev), INTENT(in) :: RHcl ! humidite relative ciel clair 59 59 INTEGER,INTENT(in) :: flag_aerosol 60 LOGICAL,INTENT(in) :: flag_bc_internal_mixture 60 61 REAL, DIMENSION(klon,klev), INTENT(in) :: pplay 61 62 REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri … … 106 107 107 108 ! 109 ! BC internal mixture 110 ! 111 INTEGER, PARAMETER :: nbclassbc = 5 ! Added by Rong Wang/OB for the 5 fractions 112 ! of BC in the soluble mode: 113 ! bc_content/0.001, 0.01, 0.02, 0.05, ! 0.1/ 114 ! for Maxwell-Garnet internal mixture 115 ! Detailed theory can be found in R. Wang Estimation of global black carbon ! direct 116 ! radiative forcing and its uncertainty constrained by observations. J. 117 ! Geophys. Res. Atmos. Added by R. Wang and OB 118 REAL :: alpha_MG_5wv(nbre_RH,las,nbclassbc) 119 120 ! 108 121 ! Proprietes optiques 109 122 ! 110 REAL :: fact_RH(nbre_RH) 111 INTEGER :: n 123 REAL :: fact_RH(nbre_RH), BC_massfra 124 INTEGER :: n, classbc 112 125 113 126 ! From here on we look at the optical parameters at 5 wavelengths: … … 117 130 DATA alpha_aers_5wv/ & 118 131 ! bc soluble 119 7.930,7.930,7.930,7.930,7.930,7.930, & 120 7.930,7.930,10.893,12.618,14.550,16.613, & 121 7.658,7.658,7.658,7.658,7.658,7.658, & 122 7.658,7.658,10.351,11.879,13.642,15.510, & 123 7.195,7.195,7.195,7.195,7.195,7.195, & 124 7.195,7.195,9.551,10.847,12.381,13.994, & 125 6.736,6.736,6.736,6.736,6.736,6.736, & 126 6.736,6.736,8.818,9.938,11.283,12.687, & 127 6.277,6.277,6.277,6.277,6.277,6.277, & 128 6.277,6.277,8.123,9.094,10.275,11.501, & 132 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 133 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 134 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 135 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 136 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501, & 129 137 ! pom soluble 130 6.676,6.676,6.676,6.676,6.710,6.934, & 131 7.141,7.569,8.034,8.529,9.456,10.511, & 132 5.109,5.109,5.109,5.109,5.189,5.535, & 133 5.960,6.852,8.008,9.712,12.897,19.676, & 134 3.718,3.718,3.718,3.718,3.779,4.042, & 135 4.364,5.052,5.956,7.314,9.896,15.688, & 136 2.849,2.849,2.849,2.849,2.897,3.107, & 137 3.365,3.916,4.649,5.760,7.900,12.863, & 138 2.229,2.229,2.229,2.229,2.268,2.437, & 139 2.645,3.095,3.692,4.608,6.391,10.633, & 138 6.676,6.676,6.676,6.676,6.710,6.934,7.141,7.569,8.034,8.529,9.456,10.511, & 139 5.109,5.109,5.109,5.109,5.189,5.535,5.960,6.852,8.008,9.712,12.897,19.676, & 140 3.718,3.718,3.718,3.718,3.779,4.042,4.364,5.052,5.956,7.314,9.896,15.688, & 141 2.849,2.849,2.849,2.849,2.897,3.107,3.365,3.916,4.649,5.760,7.900,12.863, & 142 2.229,2.229,2.229,2.229,2.268,2.437,2.645,3.095,3.692,4.608,6.391,10.633, & 140 143 ! Sulfate (Accumulation) 141 5.751,6.215,6.690,7.024,7.599,8.195, & 142 9.156,10.355,12.660,14.823,18.908,24.508, & 143 4.320,4.675,5.052,5.375,5.787,6.274, & 144 7.066,8.083,10.088,12.003,15.697,21.133, & 145 3.079,3.351,3.639,3.886,4.205,4.584, & 146 5.206,6.019,7.648,9.234,12.391,17.220, & 147 2.336,2.552,2.781,2.979,3.236,3.540, & 148 4.046,4.711,6.056,7.388,10.093,14.313, & 149 1.777,1.949,2.134,2.292,2.503,2.751, & 150 3.166,3.712,4.828,5.949,8.264,11.922, & 144 5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, & 145 4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133, & 146 3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220, & 147 2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313, & 148 1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922, & 151 149 ! Sulfate (Coarse) 152 5.751,6.215,6.690,7.024,7.599,8.195, & 153 9.156,10.355,12.660,14.823,18.908,24.508, & 154 4.320,4.675,5.052,5.375,5.787,6.274, & 155 7.066,8.083,10.088,12.003,15.697,21.133, & 156 3.079,3.351,3.639,3.886,4.205,4.584, & 157 5.206,6.019,7.648,9.234,12.391,17.220, & 158 2.336,2.552,2.781,2.979,3.236,3.540, & 159 4.046,4.711,6.056,7.388,10.093,14.313, & 160 1.777,1.949,2.134,2.292,2.503,2.751, & 161 3.166,3.712,4.828,5.949,8.264,11.922, & 150 5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, & 151 4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133, & 152 3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220, & 153 2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313, & 154 1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922, & 162 155 ! seasalt seasalt Super Coarse Soluble (SS) 163 0.218, 0.272, 0.293, 0.316, 0.343, 0.380, & 164 0.429, 0.501, 0.636, 0.755, 0.967, 1.495, & 165 0.221, 0.275, 0.297, 0.320, 0.348, 0.383, & 166 0.432, 0.509, 0.640, 0.759, 0.972, 1.510, & 167 0.224, 0.279, 0.301, 0.324, 0.352, 0.388, & 168 0.438, 0.514, 0.647, 0.768, 0.985, 1.514, & 169 0.227, 0.282, 0.303, 0.327, 0.356, 0.392, & 170 0.441, 0.518, 0.652, 0.770, 0.987, 1.529, & 171 0.230, 0.285, 0.306, 0.330, 0.359, 0.396, & 172 0.446, 0.522, 0.656, 0.777, 0.993, 1.539, & 156 0.218, 0.272, 0.293, 0.316, 0.343, 0.380, 0.429, 0.501, 0.636, 0.755, 0.967, 1.495, & 157 0.221, 0.275, 0.297, 0.320, 0.348, 0.383, 0.432, 0.509, 0.640, 0.759, 0.972, 1.510, & 158 0.224, 0.279, 0.301, 0.324, 0.352, 0.388, 0.438, 0.514, 0.647, 0.768, 0.985, 1.514, & 159 0.227, 0.282, 0.303, 0.327, 0.356, 0.392, 0.441, 0.518, 0.652, 0.770, 0.987, 1.529, & 160 0.230, 0.285, 0.306, 0.330, 0.359, 0.396, 0.446, 0.522, 0.656, 0.777, 0.993, 1.539, & 173 161 ! seasalt seasalt Coarse Soluble (CS) 174 0.578, 0.706, 0.756, 0.809, 0.876, 0.964, & 175 1.081, 1.256, 1.577, 1.858, 2.366, 3.613, & 176 0.598, 0.725, 0.779, 0.833, 0.898, 0.990, & 177 1.109, 1.290, 1.609, 1.889, 2.398, 3.682, & 178 0.619, 0.750, 0.802, 0.857, 0.927, 1.022, & 179 1.141, 1.328, 1.648, 1.939, 2.455, 3.729, & 180 0.633, 0.767, 0.820, 0.879, 0.948, 1.044, & 181 1.167, 1.353, 1.683, 1.969, 2.491, 3.785, & 182 0.648, 0.785, 0.838, 0.896, 0.967, 1.066, & 183 1.192, 1.381, 1.714, 2.006, 2.531, 3.836, & 162 0.578, 0.706, 0.756, 0.809, 0.876, 0.964, 1.081, 1.256, 1.577, 1.858, 2.366, 3.613, & 163 0.598, 0.725, 0.779, 0.833, 0.898, 0.990, 1.109, 1.290, 1.609, 1.889, 2.398, 3.682, & 164 0.619, 0.750, 0.802, 0.857, 0.927, 1.022, 1.141, 1.328, 1.648, 1.939, 2.455, 3.729, & 165 0.633, 0.767, 0.820, 0.879, 0.948, 1.044, 1.167, 1.353, 1.683, 1.969, 2.491, 3.785, & 166 0.648, 0.785, 0.838, 0.896, 0.967, 1.066, 1.192, 1.381, 1.714, 2.006, 2.531, 3.836, & 184 167 ! seasalt seasalt Accumulation Soluble (AS) 185 4.432, 5.899, 6.505, 7.166, 7.964, 7.962, & 186 9.232,11.257,14.979,18.337,24.223,37.811, & 187 3.298, 4.569, 5.110, 5.709, 6.446, 6.268, & 188 7.396, 9.246,12.787,16.113,22.197,37.136, & 189 2.340, 3.358, 3.803, 4.303, 4.928, 4.696, & 190 5.629, 7.198,10.308,13.342,19.120,34.296, & 191 1.789, 2.626, 2.999, 3.422, 3.955, 3.730, & 192 4.519, 5.864, 8.593,11.319,16.653,31.331, & 193 1.359, 2.037, 2.343, 2.693, 3.139, 2.940, & 194 3.596, 4.729, 7.076, 9.469,14.266,28.043 / 168 4.432, 5.899, 6.505, 7.166, 7.964, 7.962, 9.232,11.257,14.979,18.337,24.223,37.811, & 169 3.298, 4.569, 5.110, 5.709, 6.446, 6.268, 7.396, 9.246,12.787,16.113,22.197,37.136, & 170 2.340, 3.358, 3.803, 4.303, 4.928, 4.696, 5.629, 7.198,10.308,13.342,19.120,34.296, & 171 1.789, 2.626, 2.999, 3.422, 3.955, 3.730, 4.519, 5.864, 8.593,11.319,16.653,31.331, & 172 1.359, 2.037, 2.343, 2.693, 3.139, 2.940, 3.596, 4.729, 7.076, 9.469,14.266,28.043 / 195 173 196 174 DATA alpha_aeri_5wv/ & … … 201 179 ! pom insoluble 202 180 5.042, 3.101, 1.890, 1.294, 0.934/ 181 182 ! internal mixture 183 DATA alpha_MG_5wv/ & 184 ! bc content = 0.001 185 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 186 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 187 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 188 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 189 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501, & 190 ! bc content = 0.01 191 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 192 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 193 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 194 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 195 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501, & 196 ! bc content = 0.02 197 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 198 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 199 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 200 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 201 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501, & 202 ! bc content = 0.05 203 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 204 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 205 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 206 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 207 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501, & 208 ! bc content = 0.10 209 7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, & 210 7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, & 211 7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994, & 212 6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687, & 213 6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501 / 214 203 215 ! 204 216 ! Initialisations … … 251 263 aerosol_name(8) = id_ASSSM_phy 252 264 aerosol_name(9) = id_CIDUSTM_phy 253 aerosol_name(10) 265 aerosol_name(10)= id_CSSO4M_phy 254 266 ENDIF 255 267 … … 289 301 soluble=.TRUE. 290 302 spsol=3 291 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD303 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD 292 304 ELSEIF (aerosol_name(m).EQ.id_CSSO4M_phy) THEN 293 305 soluble=.TRUE. 294 306 spsol=4 295 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD307 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD 296 308 ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN 297 309 soluble=.TRUE. … … 325 337 IF (soluble) THEN ! For soluble aerosol 326 338 339 !--treat special case of soluble BC internal mixture 340 IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN 341 342 DO k=1, klev 343 DO i=1, klon 344 345 BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy)) 346 347 IF (BC_massfra.GE.0.10) THEN 348 classbc = 5 349 ELSEIF (BC_massfra.GE.0.05) THEN 350 classbc = 4 351 ELSEIF (BC_massfra.GE.0.02) THEN 352 classbc = 3 353 ELSEIF (BC_massfra.GE.0.01) THEN 354 classbc = 2 355 ELSE 356 classbc = 1 357 ENDIF 358 359 tau_ae5wv_int = alpha_MG_5wv(RH_num(i,k),la,classbc)+DELTA(i,k)* & 360 (alpha_MG_5wv(RH_num(i,k)+1,la,classbc) - & 361 alpha_MG_5wv(RH_num(i,k),la,classbc)) 362 tau(i,k,la,aerindex) = m_allaer(i,k,aerindex)/1.e6*zdh(i,k)*tau_ae5wv_int*fac 363 tausum(i,la,aerindex)=tausum(i,la,aerindex)+tau(i,k,la,aerindex) 364 ENDDO 365 ENDDO 366 367 !--other cases of soluble aerosols 368 ELSE 369 327 370 DO k=1, klev 328 371 DO i=1, klon … … 334 377 ENDDO 335 378 ENDDO 379 380 ENDIF 336 381 337 ELSE ! For insoluble aerosol 382 ! cases of insoluble aerosol 383 ELSE 338 384 339 385 DO k=1, klev -
LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
r2641 r2669 3 3 ! 4 4 SUBROUTINE AEROPT_6BANDS_RRTM ( & 5 pdel, m_allaer, delt,RHcl, &5 pdel, m_allaer, RHcl, & 6 6 tau_allaer, piz_allaer, & 7 7 cg_allaer, m_allaer_pi, & 8 flag_aerosol, zrho )8 flag_aerosol, flag_bc_internal_mixture, zrho ) 9 9 10 10 USE dimphy … … 28 28 ! 29 29 REAL, DIMENSION(klon,klev), INTENT(in) :: pdel 30 REAL, INTENT(in) :: delt31 30 REAL, DIMENSION(klon,klev,naero_tot), INTENT(in) :: m_allaer 32 31 REAL, DIMENSION(klon,klev,naero_tot), INTENT(in) :: m_allaer_pi 33 32 REAL, DIMENSION(klon,klev), INTENT(in) :: RHcl ! humidite relative ciel clair 34 33 INTEGER, INTENT(in) :: flag_aerosol 34 LOGICAL, INTENT(in) :: flag_bc_internal_mixture 35 35 REAL, DIMENSION(klon,klev), INTENT(in) :: zrho 36 36 ! … … 71 71 72 72 REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) :: tau_ae 73 REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) :: tau_ae_pi74 73 REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) :: piz_ae 75 74 REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) :: cg_ae 75 76 REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) :: tau_ae_pi 77 REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: piz_ae_pi 78 REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: cg_ae_pi 76 79 ! 77 80 ! Proprietes optiques … … 83 86 REAL:: piz_aers_6bands(nbre_RH,nbands_sw_rrtm,naero_soluble) !-- unit 84 87 REAL:: piz_aeri_6bands(nbands_sw_rrtm,naero_insoluble) !-- unit 85 86 INTEGER :: id 87 REAL :: tmp_var, tmp_var_pi 88 ! 89 ! BC internal mixture 90 ! 91 INTEGER, PARAMETER :: nbclassbc = 5 ! Added by Rong Wang/OB for the 5 fractions 92 ! of BC in the soluble mode: 93 ! bc_content/0.001, 0.01, 0.02, 0.05, 0.1/ 94 ! for Maxwell-Garnet internal mixture 95 ! Detailed theory can be found in R. Wang Estimation of global black carbon direct 96 ! radiative forcing and its uncertainty constrained by observations. J. 97 ! Geophys. Res. Atmos. Added by R. Wang and OB 98 REAL :: alpha_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc) 99 REAL :: cg_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc) 100 REAL :: piz_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc) 101 ! 102 INTEGER :: id, classbc, classbc_pi 103 REAL :: tmp_var, tmp_var_pi, BC_massfra, BC_massfra_pi 104 105 ! 106 REAL, PARAMETER :: tau_min = 1.e-15 107 ! REAL, PARAMETER :: tau_min = 1.e-7 88 108 89 109 !*************************************************************************** … … 268 288 0.973, 0.973, 0.972, 0.940, 0.816, 0.663 / 269 289 290 ! Added by R. Wang (July 31 2016) 291 ! properties for BC assuming Maxwell-Garnett rule and internal mixture 292 293 DATA alpha_MG_6bands/ & 294 ! bc content = 0.001 295 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 296 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 297 5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, & 298 4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, & 299 2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, & 300 0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, & 301 ! bc content = 0.01 302 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 303 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 304 5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, & 305 4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, & 306 2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, & 307 0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, & 308 ! bc content = 0.02 309 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 310 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 311 5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, & 312 4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, & 313 2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, & 314 0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, & 315 ! bc content = 0.05 316 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 317 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 318 5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, & 319 4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, & 320 2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, & 321 0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, & 322 ! bc content = 0.10 323 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 324 6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, & 325 5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, & 326 4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, & 327 2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, & 328 0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894 / 329 330 DATA cg_MG_6bands/ & 331 ! bc content = 0.001 332 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 333 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 334 0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, & 335 0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, & 336 0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, & 337 0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, & 338 ! bc content = 0.01 339 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 340 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 341 0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, & 342 0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, & 343 0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, & 344 0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, & 345 ! bc content = 0.02 346 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 347 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 348 0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, & 349 0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, & 350 0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, & 351 0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, & 352 ! bc content = 0.05 353 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 354 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 355 0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, & 356 0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, & 357 0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, & 358 0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, & 359 ! bc content = 0.10 360 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 361 0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, & 362 0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, & 363 0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, & 364 0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, & 365 0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123 / 366 367 DATA piz_MG_6bands/ & 368 ! bc content = 0.001 369 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 370 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 371 0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, & 372 0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, & 373 0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, & 374 0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, & 375 ! bc content = 0.01 376 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 377 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 378 0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, & 379 0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, & 380 0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, & 381 0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, & 382 ! bc content = 0.02 383 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 384 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 385 0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, & 386 0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, & 387 0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, & 388 0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, & 389 ! bc content = 0.05 390 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 391 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 392 0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, & 393 0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, & 394 0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, & 395 0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, & 396 ! bc content = 0.10 397 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 398 0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, & 399 0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, & 400 0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, & 401 0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, & 402 0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207 / 403 270 404 !----BEGINNING OF CALCULATIONS 271 405 … … 342 476 343 477 tau_ae(:,:,:,:)=0. 344 tau_ae_pi(:,:,:,:)=0.345 478 piz_ae(:,:,:,:)=0. 346 479 cg_ae(:,:,:,:)=0. 480 481 tau_ae_pi(:,:,:,:)=0. 482 piz_ae_pi(:,:,:,:)=0. 483 cg_ae_pi(:,:,:,:)=0. 347 484 348 485 DO m=1,nb_aer ! tau is only computed for each mass … … 357 494 soluble=.TRUE. 358 495 spsol=3 359 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD496 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD 360 497 ELSEIF (aerosol_name(m).EQ.id_CSSO4M_phy) THEN 361 498 soluble=.TRUE. 362 499 spsol=4 363 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD500 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD 364 501 ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN 365 502 soluble=.TRUE. … … 384 521 ENDIF 385 522 523 !--shortname for aerosol index 386 524 id=aerosol_name(m) 387 525 388 526 IF (soluble) THEN 389 527 390 DO k=1, klev 391 DO i=1, klon 392 tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac 393 tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac 394 395 DO inu=1,NSW 396 397 tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ & 398 delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 399 alpha_aers_6bands(RH_num(i,k),inu,spsol)) 528 !--here we treat the special case of soluble BC internal mixture with Maxwell-Garnett rule 529 IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN 530 531 DO k=1, klev 532 DO i=1, klon 533 534 tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac 535 tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac 536 537 ! Calculate the dry BC/(BC+SUL) mass ratio for all (natural+anthropogenic) aerosols 538 BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy)) 539 540 IF (BC_massfra.GE.0.10) THEN 541 classbc = 5 542 ELSEIF (BC_massfra.GE.0.05) THEN 543 classbc = 4 544 ELSEIF (BC_massfra.GE.0.02) THEN 545 classbc = 3 546 ELSEIF (BC_massfra.GE.0.01) THEN 547 classbc = 2 548 ELSE 549 classbc = 1 550 ENDIF 551 552 ! Calculate the dry BC/(BC+SUL) mass ratio for natural aerosols 553 BC_massfra_pi = m_allaer_pi(i,k,id_ASBCM_phy)/(m_allaer_pi(i,k,id_ASBCM_phy)+m_allaer_pi(i,k,id_ASSO4M_phy)) 554 555 IF (BC_massfra_pi.GE.0.10) THEN 556 classbc_pi = 5 557 ELSEIF (BC_massfra_pi.GE.0.05) THEN 558 classbc_pi = 4 559 ELSEIF (BC_massfra_pi.GE.0.02) THEN 560 classbc_pi = 3 561 ELSEIF (BC_massfra_pi.GE.0.01) THEN 562 classbc_pi = 2 563 ELSE 564 classbc_pi = 1 565 ENDIF 566 567 DO inu=1,NSW 568 569 !--all aerosols 570 tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc)+ & 571 delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc) - & 572 alpha_MG_6bands(RH_num(i,k),inu,classbc)) 400 573 401 piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + &402 delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - &403 piz_aers_6bands(RH_num(i,k),inu,spsol))574 piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc) + & 575 delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc) - & 576 piz_MG_6bands(RH_num(i,k),inu,classbc)) 404 577 405 cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + & 406 delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 407 cg_aers_6bands(RH_num(i,k),inu,spsol)) 408 409 tau_ae(i,k,id,inu) = tmp_var*tau_ae2b_int 410 tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int 411 piz_ae(i,k,id,inu) = piz_ae2b_int 412 cg_ae(i,k,id,inu) = cg_ae2b_int 413 578 cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc) + & 579 delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc) - & 580 cg_MG_6bands(RH_num(i,k),inu,classbc)) 581 582 tau_ae(i,k,id,inu) = tmp_var*tau_ae2b_int 583 piz_ae(i,k,id,inu) = piz_ae2b_int 584 cg_ae(i,k,id,inu) = cg_ae2b_int 585 586 !--natural aerosols 587 tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc_pi)+ & 588 delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) - & 589 alpha_MG_6bands(RH_num(i,k),inu,classbc_pi)) 590 591 piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc_pi) + & 592 delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) - & 593 piz_MG_6bands(RH_num(i,k),inu,classbc_pi)) 594 595 cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc_pi) + & 596 delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) - & 597 cg_MG_6bands(RH_num(i,k),inu,classbc_pi)) 598 599 tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int 600 piz_ae_pi(i,k,id,inu) = piz_ae2b_int 601 cg_ae_pi(i,k,id,inu) = cg_ae2b_int 602 603 ENDDO 414 604 ENDDO 415 605 ENDDO 416 ENDDO 606 607 !--else treat all other cases of soluble aerosols 608 ELSE 609 610 DO k=1, klev 611 DO i=1, klon 612 tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac 613 tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac 614 615 DO inu=1,NSW 616 617 tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ & 618 delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 619 alpha_aers_6bands(RH_num(i,k),inu,spsol)) 620 621 piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + & 622 delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 623 piz_aers_6bands(RH_num(i,k),inu,spsol)) 624 625 cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + & 626 delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - & 627 cg_aers_6bands(RH_num(i,k),inu,spsol)) 628 629 tau_ae(i,k,id,inu) = tmp_var*tau_ae2b_int 630 tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int 631 piz_ae(i,k,id,inu) = piz_ae2b_int 632 cg_ae(i,k,id,inu) = cg_ae2b_int 633 634 ENDDO 635 ENDDO 636 ENDDO 637 638 !--external mixture case for soluble BC 639 IF (spsol.EQ.1) THEN 640 piz_ae_pi(:,:,id,:) = piz_ae(:,:,id,:) 641 cg_ae_pi(:,:,id,:) = cg_ae(:,:,id,:) 642 ENDIF 643 644 ENDIF 417 645 418 646 ELSE ! For all aerosol insoluble components … … 443 671 DO k=1, klev 444 672 DO i=1, klon 445 !--a nthropogenicaerosol673 !--all (natural + anthropogenic) aerosol 446 674 tau_allaer(i,k,2,inu)=tau_ae(i,k,id_ASSO4M_phy,inu)+tau_ae(i,k,id_CSSO4M_phy,inu)+ & 447 675 tau_ae(i,k,id_ASBCM_phy,inu)+tau_ae(i,k,id_AIBCM_phy,inu)+ & … … 449 677 tau_ae(i,k,id_ASSSM_phy,inu)+tau_ae(i,k,id_CSSSM_phy,inu)+ & 450 678 tau_ae(i,k,id_SSSSM_phy,inu)+ tau_ae(i,k,id_CIDUSTM_phy,inu) 451 tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu), 1e-15)452 453 piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &454 tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &455 tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ &456 tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &457 tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &458 tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &459 tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &460 tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &461 tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &679 tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu),tau_min) 680 681 piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ & 682 tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ & 683 tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ & 684 tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ & 685 tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ & 686 tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ & 687 tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ & 688 tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ & 689 tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ & 462 690 tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) & 463 691 /tau_allaer(i,k,2,inu) 464 piz_allaer(i,k,2,inu)=MAX(piz_allaer(i,k,2,inu),0.01) 692 piz_allaer(i,k,2,inu)=MIN(MAX(piz_allaer(i,k,2,inu),0.01),1.0) 693 IF (tau_allaer(i,k,2,inu).LE.tau_min) piz_allaer(i,k,2,inu)=1.0 465 694 466 695 cg_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ & 467 696 tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ & 468 tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ &469 tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ &697 tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ & 698 tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ & 470 699 tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ & 471 700 tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ & 472 tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ &473 tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ &474 tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ &701 tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ & 702 tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ & 703 tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ & 475 704 tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ & 476 705 (tau_allaer(i,k,2,inu)*piz_allaer(i,k,2,inu)) 706 cg_allaer(i,k,2,inu)=MIN(MAX(cg_allaer(i,k,2,inu),0.0),1.0) 477 707 478 708 !--natural aerosol 709 !--ASBCM aerosols take _pi value because of internal mixture option 479 710 tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M_phy,inu)+tau_ae_pi(i,k,id_CSSO4M_phy,inu)+ & 480 711 tau_ae_pi(i,k,id_ASBCM_phy,inu)+tau_ae_pi(i,k,id_AIBCM_phy,inu)+ & … … 482 713 tau_ae_pi(i,k,id_ASSSM_phy,inu)+tau_ae_pi(i,k,id_CSSSM_phy,inu)+ & 483 714 tau_ae_pi(i,k,id_SSSSM_phy,inu)+ tau_ae_pi(i,k,id_CIDUSTM_phy,inu) 484 tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu), 1e-15)485 486 piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &487 tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &488 tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae (i,k,id_ASBCM_phy,inu)+&489 tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &490 tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &491 tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &492 tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &493 tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &494 tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &715 tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),tau_min) 716 717 piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ & 718 tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ & 719 tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)+ & 720 tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ & 721 tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ & 722 tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ & 723 tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ & 724 tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ & 725 tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ & 495 726 tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) & 496 727 /tau_allaer(i,k,1,inu) 497 piz_allaer(i,k,1,inu)=MAX(piz_allaer(i,k,1,inu),0.01) 498 499 cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ & 500 tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ & 501 tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ & 502 tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ & 503 tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ & 504 tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ & 505 tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ & 506 tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ & 507 tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ & 728 piz_allaer(i,k,1,inu)=MIN(MAX(piz_allaer(i,k,1,inu),0.01),1.0) 729 IF (tau_allaer(i,k,1,inu).LE.tau_min) piz_allaer(i,k,1,inu)=1.0 730 731 cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ & 732 tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ & 733 tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)*cg_ae_pi(i,k,id_ASBCM_phy,inu)+ & 734 tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ & 735 tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ & 736 tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ & 737 tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ & 738 tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ & 739 tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ & 508 740 tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ & 509 741 (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu)) 742 cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0) 510 743 511 744 ENDDO … … 513 746 ENDDO 514 747 515 !--waveband 2 and all aerosol 748 !--waveband 2 and all aerosol (third index = 2) 516 749 inu=2 517 750 DO i=1, klon -
LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosol_optic_rrtm.F90
r2542 r2669 2 2 ! 3 3 SUBROUTINE readaerosol_optic_rrtm(debut, aerosol_couple, & 4 new_aod, flag_aerosol, itap, rjourvrai, &4 new_aod, flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, & 5 5 pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & 6 6 tr_seri, mass_solu_aero, mass_solu_aero_pi, & … … 32 32 LOGICAL, INTENT(IN) :: new_aod 33 33 INTEGER, INTENT(IN) :: flag_aerosol 34 LOGICAL, INTENT(IN) :: flag_bc_internal_mixture 34 35 INTEGER, INTENT(IN) :: itap 35 36 REAL, INTENT(IN) :: rjourvrai … … 303 304 !--new aerosol properties 304 305 ! aeropt_6bands for rrtm 305 CALL aeropt_6bands_rrtm( &306 pdel, m_allaer, pdtphys, rhcl,&307 tau_aero, piz_aero, cg_aero, 308 m_allaer_pi, flag_aerosol, &309 zrho )306 CALL aeropt_6bands_rrtm( & 307 pdel, m_allaer, rhcl, & 308 tau_aero, piz_aero, cg_aero, & 309 m_allaer_pi, flag_aerosol, & 310 flag_bc_internal_mixture, zrho ) 310 311 311 312 ! aeropt_5wv only for validation and diagnostics 312 CALL aeropt_5wv_rrtm( & 313 pdel, m_allaer, & 314 pdtphys, rhcl, aerindex, & 315 flag_aerosol, pplay, t_seri, & 313 CALL aeropt_5wv_rrtm( & 314 pdel, m_allaer, & 315 rhcl, aerindex, & 316 flag_aerosol, & 317 flag_bc_internal_mixture, & 318 pplay, t_seri, & 316 319 tausum_aero, tau3d_aero ) 317 320
Note: See TracChangeset
for help on using the changeset viewer.