Changeset 2729 for LMDZ5/branches
- Timestamp:
- Dec 13, 2016, 5:13:39 PM (7 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 5 edited
-
. (modified) (1 prop)
-
DefLists/field_def_lmdz.xml (modified) (1 diff)
-
libf/phylmd/rrtm/read_rsun_rrtm.F90 (modified) (2 diffs)
-
libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 (modified) (3 diffs)
-
libf/phylmd/yamada4.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2721-2727
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/DefLists/field_def_lmdz.xml
r2720 r2729 880 880 <field id="c06_clmodis" long_name="MODIS Cloud Area Fraction" unit="1" axis_ref="pressure2" /> 881 881 <field id="c07_clmodis" long_name="MODIS Cloud Area Fraction" unit="1" axis_ref="pressure2" /> 882 <field id="c01_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 883 <field id="c02_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 884 <field id="c03_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 885 <field id="c04_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 886 <field id="c05_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 887 <field id="c06_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 888 <field id="c07_crimodis" long_name="MODIS Optical_Thickness_vs_ReffIce" unit="1" axis_ref="ReffIce" /> 889 <field id="c01_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 890 <field id="c02_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 891 <field id="c03_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 892 <field id="c04_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 893 <field id="c05_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 894 <field id="c06_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 895 <field id="c07_crlmodis" long_name="MODIS Optical_Thickness_vs_ReffLiq" unit="1" axis_ref="ReffLiq" /> 882 896 </field_group> 883 897 -
LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90
r2542 r2729 56 56 ENDIF 57 57 58 IF (size(time).NE.year_len) THEN 58 !--test if time is different from year_len but allow a mismatch of 1 day 59 IF (size(time).NE.year_len.AND.size(time).NE.year_len+1) THEN 59 60 PRINT *,'read_rsun_rrtm time <> year_len = ', size(time), year_len 60 61 CALL abort_physic('read_rsun_rrtm','time dim should be the number of days in year',1) 62 ENDIF 63 !--warning only if forcing file has 366 days but year_len has only 365 64 IF (size(time).EQ.year_len+1) THEN 65 PRINT *,'Warning read_rsun_rrtm uses a leap year rsun for a noleap year' 61 66 ENDIF 62 67 … … 91 96 92 97 !--copy 93 RSUN( :)=SSI_FRAC(:,day_cur)98 RSUN(1:NSW)=SSI_FRAC(:,day_cur) 94 99 solaire=TSI(day_cur) 95 100 -
LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
r2720 r2729 81 81 82 82 !--only root reads the data 83 IF (is_mpi_root.AND.is_omp_root) THEN83 IF (is_mpi_root.AND.is_omp_root) THEN 84 84 85 85 !--check mth_cur 86 IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN87 print *,'probleme avec le mois dans readaerosolstrat =', mth_cur88 ENDIF86 IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN 87 print *,'probleme avec le mois dans readaerosolstrat =', mth_cur 88 ENDIF 89 89 90 90 !--initialize n_lon as input data is 2D (lat-alt) only 91 n_lon = nbp_lon91 n_lon = nbp_lon 92 92 93 93 !--Starts with SW optical properties 94 94 95 CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)96 97 CALL nf95_inq_varid(ncid_in, "LEV", varid)98 CALL nf95_gw_var(ncid_in, varid, lev)99 n_lev = size(lev)100 IF (n_lev.NE.klev) THEN101 abort_message='Le nombre de niveaux n est pas egal a klev'102 CALL abort_physic(modname,abort_message,1)103 ENDIF104 105 CALL nf95_inq_varid(ncid_in, "LAT", varid)106 CALL nf95_gw_var(ncid_in, varid, latitude)107 n_lat = size(latitude)108 IF (n_lat.NE.nbp_lat) THEN109 print *, 'latitude=', n_lat, nbp_lat110 abort_message='Le nombre de lat n est pas egal a nbp_lat'111 CALL abort_physic(modname,abort_message,1)112 ENDIF113 114 CALL nf95_inq_varid(ncid_in, "TIME", varid)115 CALL nf95_gw_var(ncid_in, varid, time)116 n_month = size(time)117 IF (n_month.NE.12) THEN118 abort_message='Le nombre de month n est pas egal a 12'119 CALL abort_physic(modname,abort_message,1)120 ENDIF121 122 CALL nf95_inq_varid(ncid_in, "WAV", varid)123 CALL nf95_gw_var(ncid_in, varid, wav)124 n_wav = size(wav)125 print *, 'WAV aerosol strato=', n_wav, wav126 IF (n_wav.NE.NSW) THEN127 abort_message='Le nombre de wav n est pas egal a NSW'128 CALL abort_physic(modname,abort_message,1)129 ENDIF130 131 ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))132 ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))133 ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))134 135 ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))136 ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))137 ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))138 139 ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))140 ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))141 ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))95 CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in) 96 97 CALL nf95_inq_varid(ncid_in, "LEV", varid) 98 CALL nf95_gw_var(ncid_in, varid, lev) 99 n_lev = size(lev) 100 IF (n_lev.NE.klev) THEN 101 abort_message='Le nombre de niveaux n est pas egal a klev' 102 CALL abort_physic(modname,abort_message,1) 103 ENDIF 104 105 CALL nf95_inq_varid(ncid_in, "LAT", varid) 106 CALL nf95_gw_var(ncid_in, varid, latitude) 107 n_lat = size(latitude) 108 IF (n_lat.NE.nbp_lat) THEN 109 print *, 'latitude=', n_lat, nbp_lat 110 abort_message='Le nombre de lat n est pas egal a nbp_lat' 111 CALL abort_physic(modname,abort_message,1) 112 ENDIF 113 114 CALL nf95_inq_varid(ncid_in, "TIME", varid) 115 CALL nf95_gw_var(ncid_in, varid, time) 116 n_month = size(time) 117 IF (n_month.NE.12) THEN 118 abort_message='Le nombre de month n est pas egal a 12' 119 CALL abort_physic(modname,abort_message,1) 120 ENDIF 121 122 CALL nf95_inq_varid(ncid_in, "WAV", varid) 123 CALL nf95_gw_var(ncid_in, varid, wav) 124 n_wav = size(wav) 125 print *, 'WAV aerosol strato=', n_wav, wav 126 IF (n_wav.NE.NSW) THEN 127 abort_message='Le nombre de wav n est pas egal a NSW' 128 CALL abort_physic(modname,abort_message,1) 129 ENDIF 130 131 ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month)) 132 ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month)) 133 ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month)) 134 135 ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) 136 ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) 137 ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) 138 139 ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav)) 140 ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav)) 141 ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav)) 142 142 143 143 !--reading stratospheric aerosol tau per layer 144 CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)145 ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)146 print *,'code erreur readaerosolstrato=', ncerr, varid144 CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid) 145 ncerr = nf90_get_var(ncid_in, varid, tauaerstrat) 146 print *,'code erreur readaerosolstrato=', ncerr, varid 147 147 148 148 !--reading stratospheric aerosol omega per layer 149 CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)150 ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)151 print *,'code erreur readaerosolstrato=', ncerr, varid149 CALL nf95_inq_varid(ncid_in, "OME_SUN", varid) 150 ncerr = nf90_get_var(ncid_in, varid, pizaerstrat) 151 print *,'code erreur readaerosolstrato=', ncerr, varid 152 152 153 153 !--reading stratospheric aerosol g per layer 154 CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)155 ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)156 print *,'code erreur readaerosolstrato sw=', ncerr, varid157 158 CALL nf95_close(ncid_in)154 CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid) 155 ncerr = nf90_get_var(ncid_in, varid, cgaerstrat) 156 print *,'code erreur readaerosolstrato sw=', ncerr, varid 157 158 CALL nf95_close(ncid_in) 159 159 160 160 !--select the correct month 161 161 !--and copy into 1st longitude 162 tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)163 pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)164 cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur)162 tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur) 163 pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur) 164 cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur) 165 165 166 166 !--copy longitudes 167 DO i=2, n_lon168 tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)169 pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)170 cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:)171 ENDDO167 DO i=2, n_lon 168 tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:) 169 pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:) 170 cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:) 171 ENDDO 172 172 173 173 !---reduce to a klon_glo grid 174 DO band=1, NSW175 CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))176 CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))177 CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))178 ENDDO174 DO band=1, NSW 175 CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band)) 176 CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band)) 177 CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band)) 178 ENDDO 179 179 180 180 !--Now LW optical properties 181 181 ! 182 CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)183 184 CALL nf95_inq_varid(ncid_in, "LEV", varid)185 CALL nf95_gw_var(ncid_in, varid, lev)186 n_lev = size(lev)187 IF (n_lev.NE.klev) THEN188 abort_message='Le nombre de niveaux n est pas egal a klev'189 CALL abort_physic(modname,abort_message,1)190 ENDIF191 192 CALL nf95_inq_varid(ncid_in, "LAT", varid)193 CALL nf95_gw_var(ncid_in, varid, latitude)194 n_lat = size(latitude)195 IF (n_lat.NE.nbp_lat) THEN196 abort_message='Le nombre de lat n est pas egal a nbp_lat'197 CALL abort_physic(modname,abort_message,1)198 ENDIF199 200 CALL nf95_inq_varid(ncid_in, "TIME", varid)201 CALL nf95_gw_var(ncid_in, varid, time)202 n_month = size(time)203 IF (n_month.NE.12) THEN204 abort_message='Le nombre de month n est pas egal a 12'205 CALL abort_physic(modname,abort_message,1)206 ENDIF207 208 CALL nf95_inq_varid(ncid_in, "WAV", varid)209 CALL nf95_gw_var(ncid_in, varid, wav)210 n_wav = size(wav)211 print *, 'WAV aerosol strato=', n_wav, wav212 IF (n_wav.NE.NLW) THEN213 abort_message='Le nombre de wav n est pas egal a NLW'214 CALL abort_physic(modname,abort_message,1)215 ENDIF216 217 ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))218 ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))219 ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))182 CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in) 183 184 CALL nf95_inq_varid(ncid_in, "LEV", varid) 185 CALL nf95_gw_var(ncid_in, varid, lev) 186 n_lev = size(lev) 187 IF (n_lev.NE.klev) THEN 188 abort_message='Le nombre de niveaux n est pas egal a klev' 189 CALL abort_physic(modname,abort_message,1) 190 ENDIF 191 192 CALL nf95_inq_varid(ncid_in, "LAT", varid) 193 CALL nf95_gw_var(ncid_in, varid, latitude) 194 n_lat = size(latitude) 195 IF (n_lat.NE.nbp_lat) THEN 196 abort_message='Le nombre de lat n est pas egal a nbp_lat' 197 CALL abort_physic(modname,abort_message,1) 198 ENDIF 199 200 CALL nf95_inq_varid(ncid_in, "TIME", varid) 201 CALL nf95_gw_var(ncid_in, varid, time) 202 n_month = size(time) 203 IF (n_month.NE.12) THEN 204 abort_message='Le nombre de month n est pas egal a 12' 205 CALL abort_physic(modname,abort_message,1) 206 ENDIF 207 208 CALL nf95_inq_varid(ncid_in, "WAV", varid) 209 CALL nf95_gw_var(ncid_in, varid, wav) 210 n_wav = size(wav) 211 print *, 'WAV aerosol strato=', n_wav, wav 212 IF (n_wav.NE.NLW) THEN 213 abort_message='Le nombre de wav n est pas egal a NLW' 214 CALL abort_physic(modname,abort_message,1) 215 ENDIF 216 217 ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month)) 218 ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) 219 ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav)) 220 220 221 221 !--reading stratospheric aerosol lw tau per layer 222 CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)223 ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)224 print *,'code erreur readaerosolstrato lw=', ncerr, varid225 226 CALL nf95_close(ncid_in)222 CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid) 223 ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat) 224 print *,'code erreur readaerosolstrato lw=', ncerr, varid 225 226 CALL nf95_close(ncid_in) 227 227 228 228 !--select the correct month 229 229 !--and copy into 1st longitude 230 taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)230 taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur) 231 231 !--copy longitudes 232 DO i=2, n_lon233 taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)234 ENDDO232 DO i=2, n_lon 233 taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:) 234 ENDDO 235 235 236 236 !---reduce to a klon_glo grid 237 DO band=1, NLW 238 CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band)) 239 ENDDO 240 241 ENDIF !--is_mpi_root and is_omp_root 237 DO band=1, NLW 238 CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band)) 239 ENDDO 240 241 ELSE !--proc other than mpi_root and omp_root 242 !--dummy allocation needed for debug mode 243 244 ALLOCATE(tauaerstrat_mois_glo(1,1,1)) 245 ALLOCATE(pizaerstrat_mois_glo(1,1,1)) 246 ALLOCATE(cgaerstrat_mois_glo(1,1,1)) 247 ALLOCATE(taulwaerstrat_mois_glo(1,1,1)) 248 249 ENDIF !--is_mpi_root and is_omp_root 242 250 243 251 !$OMP BARRIER 244 252 245 253 !--keep memory of previous month 246 mth_pre=mth_cur254 mth_pre=mth_cur 247 255 248 256 !--scatter on all proc 249 CALL scatter(tauaerstrat_mois_glo,tau_aer_strat) 250 CALL scatter(pizaerstrat_mois_glo,piz_aer_strat) 251 CALL scatter(cgaerstrat_mois_glo,cg_aer_strat) 252 CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat) 253 254 IF (is_mpi_root.AND.is_omp_root) THEN 255 ! 256 DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat) 257 DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois) 258 DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo) 259 260 DEALLOCATE(taulwaerstrat,taulwaerstrat_mois,taulwaerstrat_mois_glo) 261 ! 262 ENDIF !--is_mpi_root and is_omp_root 257 CALL scatter(tauaerstrat_mois_glo,tau_aer_strat) 258 CALL scatter(pizaerstrat_mois_glo,piz_aer_strat) 259 CALL scatter(cgaerstrat_mois_glo,cg_aer_strat) 260 CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat) 261 262 IF (is_mpi_root.AND.is_omp_root) THEN 263 ! 264 DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat) 265 DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois) 266 DEALLOCATE(taulwaerstrat,taulwaerstrat_mois) 267 ! 268 ENDIF !--is_mpi_root and is_omp_root 269 270 DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo) 271 DEALLOCATE(taulwaerstrat_mois_glo) 263 272 264 273 !$OMP BARRIER … … 282 291 !--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones 283 292 DO band=1, NSW 284 WHERE (stratomask.GT.0.999999)293 WHERE (stratomask.GT.0.999999) 285 294 !--anthropogenic aerosols bands 1 to NSW 286 cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &287 cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &288 MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &289 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )290 piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &291 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &292 MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )293 tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)295 cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & 296 cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & 297 MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & 298 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) 299 piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & 300 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & 301 MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 ) 302 tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band) 294 303 !--natural aerosols bands 1 to NSW 295 cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &296 cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &297 MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &298 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )299 piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &300 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &301 MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )302 tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)304 cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & 305 cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & 306 MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & 307 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) 308 piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & 309 piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & 310 MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 ) 311 tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band) 303 312 !--no stratospheric aerosol in index 1 for these tests 304 ! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band)305 ! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band)306 ! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band)313 ! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band) 314 ! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band) 315 ! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) 307 316 ENDWHERE 308 317 ENDDO … … 322 331 323 332 DO band=1, NLW 324 WHERE (stratomask.GT.0.999999)325 tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)326 tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)333 WHERE (stratomask.GT.0.999999) 334 tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band) 335 tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band) 327 336 !--no stratospheric aerosols in index 1 for these tests 328 337 ! tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) 329 ENDWHERE338 ENDWHERE 330 339 ENDDO 331 340 -
LMDZ5/branches/testing/libf/phylmd/yamada4.F90
r2594 r2729 2 2 3 3 SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, & 4 cd, q2, km, kn, kq, ustar, iflag_pbl)4 cd, tke, km, kn, kq, ustar, iflag_pbl) 5 5 6 6 USE dimphy … … 56 56 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l 57 57 ! -> the model can run with longer time-steps. 58 ! 2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr) 59 ! On met tke (=q2/2) en entrée plutôt que q2 60 ! On corrige l'update de la tke 58 61 ! 59 62 ! Inpout/Output : 60 63 !============== 61 ! q2 : $q^2$au bas de chaque couche64 ! tke : tke au bas de chaque couche 62 65 ! (en entree : la valeur au debut du pas de temps) 63 66 ! (en sortie : la valeur a la fin du pas de temps) … … 90 93 REAL teta(klon, klev) 91 94 REAL cd(klon) 92 REAL q2(klon, klev+1)95 REAL tke(klon, klev+1) 93 96 REAL unsdz(klon, klev) 94 97 REAL unsdzdec(klon, klev+1) … … 104 107 INCLUDE "clesphys.h" 105 108 106 109 REAL q2(klon, klev+1) 107 110 REAL kmpre(klon, klev+1), tmp2, qpre 108 111 REAL mpre(klon, klev+1) … … 127 130 REAL zz(klon, klev+1) 128 131 INTEGER iter 129 REAL ric,rifc, b1, kap130 SAVE ric, rifc, b1, kap131 DATA ric, rifc, b1, kap/0.195, 0.191,16.6, 0.4/132 !$OMP THREADPRIVATE(ric,rifc,b1,kap)133 REAL seuilsm, seuilalpha132 REAL,SAVE :: ric0,ric,rifc, b1, kap 133 !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap) 134 DATA b1, kap/16.6, 0.4/ 135 REAL,SAVE :: seuilsm, seuilalpha 136 !$OMP THREADPRIVATE(seuilsm, seuilalpha) 134 137 REAL,SAVE :: lmixmin 135 138 !$OMP THREADPRIVATE(lmixmin) 139 LOGICAL, SAVE :: new_yamada4 140 !$OMP THREADPRIVATE(new_yamada4) 141 REAL, SAVE :: yun,ydeux 142 !$OMP THREADPRIVATE(yun,ydeux) 136 143 REAL frif, falpha, fsm 137 144 REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), & … … 141 148 142 149 150 143 151 ! Fonctions utilis??es 144 152 !-------------------- … … 150 158 151 159 IF (firstcall) THEN 160 ! Seuil dans le code de turbulence 161 new_yamada4=.false. 162 CALL getin_p('new_yamada4',new_yamada4) 163 ! Pour garantir la continuite dans la mise au point de CMIP6. 164 ! Eliminer l'option new_yamada4 des le printemps 2016 165 IF (new_yamada4) THEN 166 ric=0.143 ! qui donne des valeurs proches des seuils proposes 167 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83) 168 ! sm=1.1213 (1.12 dans Y83) 169 CALL getin_p('yamada4_ric',ric) 170 ric0=0.19489 ! ric=0.195 originalement, mais produisait sm<0 171 ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif 172 rifc=frif(ric) 173 seuilsm=fsm(frif(ric)) 174 seuilalpha=falpha(frif(ric)) 175 yun=1. 176 ydeux=2. 177 ! yun=2. 178 ! ydeux=1. 179 ELSE 180 !! DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/ 181 ric=0.195 182 rifc=0.191 183 seuilalpha=1.12 184 seuilsm=0.085 185 yun=2. 186 ydeux=1. 187 ENDIF 188 PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha 152 189 firstcall = .FALSE. 153 190 lmixmin=1. … … 169 206 END IF 170 207 171 ! Seuil dans le code de turbulence172 seuilalpha=1.12173 seuilsm=0.085174 208 175 209 nlay = klev … … 231 265 rif(ig, k) = rifc 232 266 END IF 233 267 if (new_yamada4) then 268 alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha) 269 sm(ig, k) = max(fsm(rif(ig,k)),seuilsm) 270 else 234 271 IF (rif(ig,k)<0.16) THEN 235 272 alpha(ig, k) = falpha(rif(ig,k)) … … 239 276 sm(ig, k) = seuilsm 240 277 END IF 278 279 end if 241 280 zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k) 242 281 END DO … … 244 283 245 284 285 286 287 288 !======================================================================= 289 ! DIFFERENT TYPES DE SCHEMA de YAMADA 290 !======================================================================= 291 292 ! On commence par calculé q2 à partir de la tke 293 294 IF (new_yamada4) THEN 295 DO k=1,klev+1 296 q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux 297 ENDDO 298 ELSE 299 DO k=1,klev+1 300 q2(1:ngrid,k)=tke(1:ngrid,k) 301 ENDDO 302 ENDIF 246 303 247 304 ! ==================================================================== … … 252 309 CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l) 253 310 254 255 256 !=======================================================================257 ! DIFFERENT TYPES DE SCHEMA de YAMADA258 !=======================================================================259 311 260 312 !-------------- … … 350 402 DO k = 2, klev - 1 351 403 km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) 352 q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) 404 q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) 405 ! q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) 353 406 q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) 354 q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) 407 q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1)) 408 ! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) 355 409 q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) 356 410 END DO … … 517 571 lyam(1:ngrid, 2:klev) 518 572 END IF 573 574 575 576 !============================================================================ 577 ! Mise à jour de la tke 578 !============================================================================ 579 580 IF (new_yamada4) THEN 581 DO k=1,klev+1 582 tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux 583 ENDDO 584 ELSE 585 DO k=1,klev+1 586 tke(1:ngrid,k)=q2(1:ngrid,k) 587 ENDDO 588 ENDIF 589 519 590 520 591 !============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
