Index: LMDZ5/branches/testing/DefLists/field_def_lmdz.xml =================================================================== --- LMDZ5/branches/testing/DefLists/field_def_lmdz.xml (revision 2720) +++ LMDZ5/branches/testing/DefLists/field_def_lmdz.xml (revision 2729) @@ -880,4 +880,18 @@ + + + + + + + + + + + + + + Index: LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90 =================================================================== --- LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90 (revision 2720) +++ LMDZ5/branches/testing/libf/phylmd/rrtm/read_rsun_rrtm.F90 (revision 2729) @@ -56,7 +56,12 @@ ENDIF - IF (size(time).NE.year_len) THEN +!--test if time is different from year_len but allow a mismatch of 1 day + IF (size(time).NE.year_len.AND.size(time).NE.year_len+1) THEN PRINT *,'read_rsun_rrtm time <> year_len = ', size(time), year_len CALL abort_physic('read_rsun_rrtm','time dim should be the number of days in year',1) + ENDIF +!--warning only if forcing file has 366 days but year_len has only 365 + IF (size(time).EQ.year_len+1) THEN + PRINT *,'Warning read_rsun_rrtm uses a leap year rsun for a noleap year' ENDIF @@ -91,5 +96,5 @@ !--copy - RSUN(:)=SSI_FRAC(:,day_cur) + RSUN(1:NSW)=SSI_FRAC(:,day_cur) solaire=TSI(day_cur) Index: LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 =================================================================== --- LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 (revision 2720) +++ LMDZ5/branches/testing/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 (revision 2729) @@ -81,184 +81,193 @@ !--only root reads the data - IF (is_mpi_root.AND.is_omp_root) THEN + IF (is_mpi_root.AND.is_omp_root) THEN !--check mth_cur - IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN - print *,'probleme avec le mois dans readaerosolstrat =', mth_cur - ENDIF + IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN + print *,'probleme avec le mois dans readaerosolstrat =', mth_cur + ENDIF !--initialize n_lon as input data is 2D (lat-alt) only - n_lon = nbp_lon + n_lon = nbp_lon !--Starts with SW optical properties - CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in) - - CALL nf95_inq_varid(ncid_in, "LEV", varid) - CALL nf95_gw_var(ncid_in, varid, lev) - n_lev = size(lev) - IF (n_lev.NE.klev) THEN - abort_message='Le nombre de niveaux n est pas egal a klev' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "LAT", varid) - CALL nf95_gw_var(ncid_in, varid, latitude) - n_lat = size(latitude) - IF (n_lat.NE.nbp_lat) THEN - print *, 'latitude=', n_lat, nbp_lat - abort_message='Le nombre de lat n est pas egal a nbp_lat' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "TIME", varid) - CALL nf95_gw_var(ncid_in, varid, time) - n_month = size(time) - IF (n_month.NE.12) THEN - abort_message='Le nombre de month n est pas egal a 12' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "WAV", varid) - CALL nf95_gw_var(ncid_in, varid, wav) - n_wav = size(wav) - print *, 'WAV aerosol strato=', n_wav, wav - IF (n_wav.NE.NSW) THEN - abort_message='Le nombre de wav n est pas egal a NSW' - CALL abort_physic(modname,abort_message,1) - ENDIF - - ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month)) - ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month)) - ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month)) - - ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) - ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) - ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) - - ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav)) - ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav)) - ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav)) + CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in) + + CALL nf95_inq_varid(ncid_in, "LEV", varid) + CALL nf95_gw_var(ncid_in, varid, lev) + n_lev = size(lev) + IF (n_lev.NE.klev) THEN + abort_message='Le nombre de niveaux n est pas egal a klev' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "LAT", varid) + CALL nf95_gw_var(ncid_in, varid, latitude) + n_lat = size(latitude) + IF (n_lat.NE.nbp_lat) THEN + print *, 'latitude=', n_lat, nbp_lat + abort_message='Le nombre de lat n est pas egal a nbp_lat' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "TIME", varid) + CALL nf95_gw_var(ncid_in, varid, time) + n_month = size(time) + IF (n_month.NE.12) THEN + abort_message='Le nombre de month n est pas egal a 12' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "WAV", varid) + CALL nf95_gw_var(ncid_in, varid, wav) + n_wav = size(wav) + print *, 'WAV aerosol strato=', n_wav, wav + IF (n_wav.NE.NSW) THEN + abort_message='Le nombre de wav n est pas egal a NSW' + CALL abort_physic(modname,abort_message,1) + ENDIF + + ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month)) + ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month)) + ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month)) + + ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) + ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) + ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) + + ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav)) + ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav)) + ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav)) !--reading stratospheric aerosol tau per layer - CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid) - ncerr = nf90_get_var(ncid_in, varid, tauaerstrat) - print *,'code erreur readaerosolstrato=', ncerr, varid + CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid) + ncerr = nf90_get_var(ncid_in, varid, tauaerstrat) + print *,'code erreur readaerosolstrato=', ncerr, varid !--reading stratospheric aerosol omega per layer - CALL nf95_inq_varid(ncid_in, "OME_SUN", varid) - ncerr = nf90_get_var(ncid_in, varid, pizaerstrat) - print *,'code erreur readaerosolstrato=', ncerr, varid + CALL nf95_inq_varid(ncid_in, "OME_SUN", varid) + ncerr = nf90_get_var(ncid_in, varid, pizaerstrat) + print *,'code erreur readaerosolstrato=', ncerr, varid !--reading stratospheric aerosol g per layer - CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid) - ncerr = nf90_get_var(ncid_in, varid, cgaerstrat) - print *,'code erreur readaerosolstrato sw=', ncerr, varid - - CALL nf95_close(ncid_in) + CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid) + ncerr = nf90_get_var(ncid_in, varid, cgaerstrat) + print *,'code erreur readaerosolstrato sw=', ncerr, varid + + CALL nf95_close(ncid_in) !--select the correct month !--and copy into 1st longitude - tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur) - pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur) - cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur) + tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur) + pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur) + cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur) !--copy longitudes - DO i=2, n_lon - tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:) - pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:) - cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:) - ENDDO + DO i=2, n_lon + tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:) + pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:) + cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:) + ENDDO !---reduce to a klon_glo grid - DO band=1, NSW - CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band)) - CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band)) - CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band)) - ENDDO + DO band=1, NSW + CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band)) + CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band)) + CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band)) + ENDDO !--Now LW optical properties ! - CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in) - - CALL nf95_inq_varid(ncid_in, "LEV", varid) - CALL nf95_gw_var(ncid_in, varid, lev) - n_lev = size(lev) - IF (n_lev.NE.klev) THEN - abort_message='Le nombre de niveaux n est pas egal a klev' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "LAT", varid) - CALL nf95_gw_var(ncid_in, varid, latitude) - n_lat = size(latitude) - IF (n_lat.NE.nbp_lat) THEN - abort_message='Le nombre de lat n est pas egal a nbp_lat' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "TIME", varid) - CALL nf95_gw_var(ncid_in, varid, time) - n_month = size(time) - IF (n_month.NE.12) THEN - abort_message='Le nombre de month n est pas egal a 12' - CALL abort_physic(modname,abort_message,1) - ENDIF - - CALL nf95_inq_varid(ncid_in, "WAV", varid) - CALL nf95_gw_var(ncid_in, varid, wav) - n_wav = size(wav) - print *, 'WAV aerosol strato=', n_wav, wav - IF (n_wav.NE.NLW) THEN - abort_message='Le nombre de wav n est pas egal a NLW' - CALL abort_physic(modname,abort_message,1) - ENDIF - - ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month)) - ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) - ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav)) + CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in) + + CALL nf95_inq_varid(ncid_in, "LEV", varid) + CALL nf95_gw_var(ncid_in, varid, lev) + n_lev = size(lev) + IF (n_lev.NE.klev) THEN + abort_message='Le nombre de niveaux n est pas egal a klev' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "LAT", varid) + CALL nf95_gw_var(ncid_in, varid, latitude) + n_lat = size(latitude) + IF (n_lat.NE.nbp_lat) THEN + abort_message='Le nombre de lat n est pas egal a nbp_lat' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "TIME", varid) + CALL nf95_gw_var(ncid_in, varid, time) + n_month = size(time) + IF (n_month.NE.12) THEN + abort_message='Le nombre de month n est pas egal a 12' + CALL abort_physic(modname,abort_message,1) + ENDIF + + CALL nf95_inq_varid(ncid_in, "WAV", varid) + CALL nf95_gw_var(ncid_in, varid, wav) + n_wav = size(wav) + print *, 'WAV aerosol strato=', n_wav, wav + IF (n_wav.NE.NLW) THEN + abort_message='Le nombre de wav n est pas egal a NLW' + CALL abort_physic(modname,abort_message,1) + ENDIF + + ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month)) + ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav)) + ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav)) !--reading stratospheric aerosol lw tau per layer - CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid) - ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat) - print *,'code erreur readaerosolstrato lw=', ncerr, varid - - CALL nf95_close(ncid_in) + CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid) + ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat) + print *,'code erreur readaerosolstrato lw=', ncerr, varid + + CALL nf95_close(ncid_in) !--select the correct month !--and copy into 1st longitude - taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur) + taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur) !--copy longitudes - DO i=2, n_lon - taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:) - ENDDO + DO i=2, n_lon + taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:) + ENDDO !---reduce to a klon_glo grid - DO band=1, NLW - CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band)) - ENDDO - - ENDIF !--is_mpi_root and is_omp_root + DO band=1, NLW + CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band)) + ENDDO + + ELSE !--proc other than mpi_root and omp_root + !--dummy allocation needed for debug mode + + ALLOCATE(tauaerstrat_mois_glo(1,1,1)) + ALLOCATE(pizaerstrat_mois_glo(1,1,1)) + ALLOCATE(cgaerstrat_mois_glo(1,1,1)) + ALLOCATE(taulwaerstrat_mois_glo(1,1,1)) + + ENDIF !--is_mpi_root and is_omp_root !$OMP BARRIER !--keep memory of previous month - mth_pre=mth_cur + mth_pre=mth_cur !--scatter on all proc - CALL scatter(tauaerstrat_mois_glo,tau_aer_strat) - CALL scatter(pizaerstrat_mois_glo,piz_aer_strat) - CALL scatter(cgaerstrat_mois_glo,cg_aer_strat) - CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat) - - IF (is_mpi_root.AND.is_omp_root) THEN -! - DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat) - DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois) - DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo) - - DEALLOCATE(taulwaerstrat,taulwaerstrat_mois,taulwaerstrat_mois_glo) -! - ENDIF !--is_mpi_root and is_omp_root + CALL scatter(tauaerstrat_mois_glo,tau_aer_strat) + CALL scatter(pizaerstrat_mois_glo,piz_aer_strat) + CALL scatter(cgaerstrat_mois_glo,cg_aer_strat) + CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat) + + IF (is_mpi_root.AND.is_omp_root) THEN +! + DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat) + DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois) + DEALLOCATE(taulwaerstrat,taulwaerstrat_mois) +! + ENDIF !--is_mpi_root and is_omp_root + + DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo) + DEALLOCATE(taulwaerstrat_mois_glo) !$OMP BARRIER @@ -282,27 +291,27 @@ !--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones DO band=1, NSW - WHERE (stratomask.GT.0.999999) + WHERE (stratomask.GT.0.999999) !--anthropogenic aerosols bands 1 to NSW - cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & - cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & - MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & - piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) - piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & - piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & - MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 ) - tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band) + cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & + cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & + MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & + piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) + piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + & + piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & + MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 ) + tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band) !--natural aerosols bands 1 to NSW - cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & - cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & - MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & - piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) - piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & - piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & - MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 ) - tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band) + cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & + cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & + MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & + piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 ) + piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + & + piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / & + MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 ) + tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band) !--no stratospheric aerosol in index 1 for these tests -! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band) -! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band) -! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) +! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band) +! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band) +! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) ENDWHERE ENDDO @@ -322,10 +331,10 @@ DO band=1, NLW - WHERE (stratomask.GT.0.999999) - tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band) - tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band) + WHERE (stratomask.GT.0.999999) + tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band) + tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band) !--no stratospheric aerosols in index 1 for these tests ! tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) - ENDWHERE + ENDWHERE ENDDO Index: LMDZ5/branches/testing/libf/phylmd/yamada4.F90 =================================================================== --- LMDZ5/branches/testing/libf/phylmd/yamada4.F90 (revision 2720) +++ LMDZ5/branches/testing/libf/phylmd/yamada4.F90 (revision 2729) @@ -2,5 +2,5 @@ SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, & - cd, q2, km, kn, kq, ustar, iflag_pbl) + cd, tke, km, kn, kq, ustar, iflag_pbl) USE dimphy @@ -56,8 +56,11 @@ ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l ! -> the model can run with longer time-steps. + ! 2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr) + ! On met tke (=q2/2) en entrée plutôt que q2 + ! On corrige l'update de la tke ! ! Inpout/Output : !============== - ! q2 : $q^2$ au bas de chaque couche + ! tke : tke au bas de chaque couche ! (en entree : la valeur au debut du pas de temps) ! (en sortie : la valeur a la fin du pas de temps) @@ -90,5 +93,5 @@ REAL teta(klon, klev) REAL cd(klon) - REAL q2(klon, klev+1) + REAL tke(klon, klev+1) REAL unsdz(klon, klev) REAL unsdzdec(klon, klev+1) @@ -104,5 +107,5 @@ INCLUDE "clesphys.h" - + REAL q2(klon, klev+1) REAL kmpre(klon, klev+1), tmp2, qpre REAL mpre(klon, klev+1) @@ -127,11 +130,15 @@ REAL zz(klon, klev+1) INTEGER iter - REAL ric, rifc, b1, kap - SAVE ric, rifc, b1, kap - DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/ - !$OMP THREADPRIVATE(ric,rifc,b1,kap) - REAL seuilsm, seuilalpha + REAL,SAVE :: ric0,ric,rifc, b1, kap + !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap) + DATA b1, kap/16.6, 0.4/ + REAL,SAVE :: seuilsm, seuilalpha + !$OMP THREADPRIVATE(seuilsm, seuilalpha) REAL,SAVE :: lmixmin !$OMP THREADPRIVATE(lmixmin) + LOGICAL, SAVE :: new_yamada4 + !$OMP THREADPRIVATE(new_yamada4) + REAL, SAVE :: yun,ydeux + !$OMP THREADPRIVATE(yun,ydeux) REAL frif, falpha, fsm REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), & @@ -141,4 +148,5 @@ + ! Fonctions utilis??es !-------------------- @@ -150,4 +158,33 @@ IF (firstcall) THEN +! Seuil dans le code de turbulence + new_yamada4=.false. + CALL getin_p('new_yamada4',new_yamada4) +! Pour garantir la continuite dans la mise au point de CMIP6. +! Eliminer l'option new_yamada4 des le printemps 2016 + IF (new_yamada4) THEN + ric=0.143 ! qui donne des valeurs proches des seuils proposes + ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83) + ! sm=1.1213 (1.12 dans Y83) + CALL getin_p('yamada4_ric',ric) + ric0=0.19489 ! ric=0.195 originalement, mais produisait sm<0 + ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif + rifc=frif(ric) + seuilsm=fsm(frif(ric)) + seuilalpha=falpha(frif(ric)) + yun=1. + ydeux=2. +! yun=2. +! ydeux=1. + ELSE + !! DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/ + ric=0.195 + rifc=0.191 + seuilalpha=1.12 + seuilsm=0.085 + yun=2. + ydeux=1. + ENDIF + PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha firstcall = .FALSE. lmixmin=1. @@ -169,7 +206,4 @@ END IF -! Seuil dans le code de turbulence - seuilalpha=1.12 - seuilsm=0.085 nlay = klev @@ -231,5 +265,8 @@ rif(ig, k) = rifc END IF - +if (new_yamada4) then + alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha) + sm(ig, k) = max(fsm(rif(ig,k)),seuilsm) +else IF (rif(ig,k)<0.16) THEN alpha(ig, k) = falpha(rif(ig,k)) @@ -239,4 +276,6 @@ sm(ig, k) = seuilsm END IF + +end if zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k) END DO @@ -244,4 +283,22 @@ + + + + !======================================================================= + ! DIFFERENT TYPES DE SCHEMA de YAMADA + !======================================================================= + + ! On commence par calculé q2 à partir de la tke + + IF (new_yamada4) THEN + DO k=1,klev+1 + q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux + ENDDO + ELSE + DO k=1,klev+1 + q2(1:ngrid,k)=tke(1:ngrid,k) + ENDDO + ENDIF ! ==================================================================== @@ -252,9 +309,4 @@ CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l) - - - !======================================================================= - ! DIFFERENT TYPES DE SCHEMA de YAMADA - !======================================================================= !-------------- @@ -350,7 +402,9 @@ DO k = 2, klev - 1 km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) - q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) + q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) +! q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) - q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) + q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1)) +! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) END DO @@ -517,4 +571,21 @@ lyam(1:ngrid, 2:klev) END IF + + + +!============================================================================ +! Mise à jour de la tke +!============================================================================ + + IF (new_yamada4) THEN + DO k=1,klev+1 + tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux + ENDDO + ELSE + DO k=1,klev+1 + tke(1:ngrid,k)=q2(1:ngrid,k) + ENDDO + ENDIF + !============================================================================