Index: LMDZ6/trunk/libf/phylmd/comsoil.h =================================================================== --- LMDZ6/trunk/libf/phylmd/comsoil.h (revision 3967) +++ LMDZ6/trunk/libf/phylmd/comsoil.h (revision 3974) @@ -4,6 +4,6 @@ common /comsoil/inertie_sol,inertie_sno,inertie_sic,inertie_lic, & - & iflag_sic + & iflag_sic,iflag_inertie real inertie_sol,inertie_sno,inertie_sic,inertie_lic - integer iflag_sic + integer iflag_sic,iflag_inertie !$OMP THREADPRIVATE(/comsoil/) Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 3967) +++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 3974) @@ -183,5 +183,5 @@ REAL,SAVE :: exposant_glace_omp REAL,SAVE :: rei_min_omp, rei_max_omp - INTEGER,SAVE :: iflag_sic_omp + INTEGER,SAVE :: iflag_sic_omp, iflag_inertie_omp REAL,SAVE :: inertie_sol_omp,inertie_sno_omp,inertie_sic_omp REAL,SAVE :: inertie_lic_omp @@ -1325,4 +1325,12 @@ CALL getin('iflag_sic',iflag_sic_omp) ! + !Config Key = iflag_inertie + !Config Desc = + !Config Def = 0 + !Config Help = + ! + iflag_inertie_omp = 0 + CALL getin('iflag_inertie',iflag_inertie_omp) + ! !Config Key = inertie_sic !Config Desc = @@ -2395,4 +2403,5 @@ albsno0 = albsno0_omp iflag_sic = iflag_sic_omp + iflag_inertie = iflag_inertie_omp inertie_sol = inertie_sol_omp inertie_sic = inertie_sic_omp @@ -2881,4 +2890,5 @@ WRITE(lunout,*) ' albsno0 = ', albsno0 WRITE(lunout,*) ' iflag_sic = ', iflag_sic + WRITE(lunout,*) ' iflag_inertie = ', iflag_inertie WRITE(lunout,*) ' inertie_sol = ', inertie_sol WRITE(lunout,*) ' inertie_sic = ', inertie_sic Index: LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90 (revision 3967) +++ LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90 (revision 3974) @@ -180,4 +180,5 @@ ! USE dimphy + USE geometry_mod, ONLY: longitude,latitude USE calcul_fluxs_mod USE surface_data, ONLY : calice, calsno @@ -260,5 +261,6 @@ IF (soil_model) THEN ! update tsoil and calculate soilcap and soilflux - CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, tsoil,soilcap, soilflux) + CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, & + & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil,soilcap, soilflux) cal(1:knon) = RCPD / soilcap(1:knon) radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) Index: LMDZ6/trunk/libf/phylmd/soil.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/soil.F90 (revision 3967) +++ LMDZ6/trunk/libf/phylmd/soil.F90 (revision 3974) @@ -2,6 +2,6 @@ ! $Header$ ! -SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, & - ptsoil, pcapcal, pfluxgrd) +SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, & + lon, lat, ptsoil, pcapcal, pfluxgrd) USE dimphy @@ -21,4 +21,8 @@ ! the surface conduction flux pcapcal ! +! Update: 2021/07 : soil thermal inertia, formerly a constant value, +! ------ can also be now a function of soil moisture (F Cheruy's idea) +! depending on iflag_inertie, read from physiq.def via conf_phys_m.F90 +! ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy) ! ! Method: Implicit time integration @@ -48,4 +52,7 @@ ! snow(klon) snow ! ptsrf(klon) surface temperature at time-step t (K) +! qsol(klon) soil moisture (kg/m2 or mm) +! lon(klon) longitude in radian +! lat(klon) latitude in radian ! ptsoil(klon,nsoilmx) temperature inside the ground (K) ! pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) @@ -60,8 +67,11 @@ ! --------- REAL, INTENT(IN) :: ptimestep - INTEGER, INTENT(IN) :: indice, knon + INTEGER, INTENT(IN) :: indice, knon !, knindex REAL, DIMENSION(klon), INTENT(IN) :: snow REAL, DIMENSION(klon), INTENT(IN) :: ptsrf - + REAL, DIMENSION(klon), INTENT(IN) :: qsol + REAL, DIMENSION(klon), INTENT(IN) :: lon + REAL, DIMENSION(klon), INTENT(IN) :: lat + REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil REAL, DIMENSION(klon), INTENT(OUT) :: pcapcal @@ -182,7 +192,52 @@ ! knon, knindex, ztherm_i) ELSE IF (indice == is_ter) THEN + ! + ! La relation entre l'inertie thermique du sol et qsol change d'apres + ! iflag_inertie, defini dans physiq.def, et appele via comsoil.h + ! DO ig = 1, knon - ztherm_i(ig) = inertie_sol + ! iflag_inertie=0 correspond au cas inertie=constant, comme avant + IF (iflag_inertie==0) THEN + ztherm_i(ig) = inertie_sol + ELSE IF (iflag_inertie == 1) THEN + ! I = a_qsol * qsol + b modele lineaire deduit d'une + ! regression lineaire I = a_mrsos * mrsos + b obtenue sur + ! sorties MO d'une simulation LMDZOR(CMIP6) sur l'annee 2000 + ! sur tous les points avec frac_snow=0 + ! Difference entre qsol et mrsos prise en compte par un + ! facteur d'echelle sur le coefficient directeur de regression: + ! fact = 35./150. = mrsos_max/qsol_max + ! et a_qsol = a_mrsos * fact (car a = dI/dHumidite) + ztherm_i(ig) = 30.0 *35.0/150.0 *qsol(ig) +770.0 + ! AS : pour qsol entre 0 - 150, on a I entre 770 - 1820 + ELSE IF (iflag_inertie == 2) THEN + ! deux regressions lineaires, sur les memes sorties, + ! distinguant le type de sol : sable ou autre (limons/argile) + ! Implementation simple : regression type "sable" seulement pour + ! Sahara, defini par une "boite" lat/lon (NB : en radians !! ) + IF (lon(ig)>-0.35 .AND. lon(ig)<0.70 .AND. lat(ig)>0.17 .AND. lat(ig)<0.52) THEN + ! Valeurs theoriquement entre 728 et 2373 ; qsol valeurs basses + ztherm_i(ig) = 47. *35.0/150.0 *qsol(ig) +728. ! boite type "sable" pour Sahara + ELSE + ! Valeurs theoriquement entre 550 et 1940 ; qsol valeurs moyennes et hautes + ztherm_i(ig) = 41. *35.0/150.0 *qsol(ig) +505. + ENDIF + ELSE IF (iflag_inertie == 3) THEN + ! AS : idee a tester : + ! si la relation doit etre une droite, + ! definissons-la en fonction des valeurs min et max de qsol (0:150), + ! et de l'inertie (900 : 2000 ou 2400 ; choix ici: 2000) + ! I = I_min + qsol * (I_max - I_min)/(qsol_max - qsol_min) + ztherm_i(ig) = 900. + qsol(ig) * (2000. - 900.)/150. + ELSE + WRITE (lunout,*) "Le choix iflag_inertie = ",iflag_inertie," n'est pas defini. Veuillez choisir un entier entre 0 et 3" + ENDIF + ! + ! Fin de l'introduction de la relation entre l'inertie thermique du sol et qsol + !------------------------------------------- + !AS : donc le moindre flocon de neige sur un point de grid + ! fait que l'inertie du point passe a la valeur pour neige ! IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno + ENDDO ! CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', & Index: LMDZ6/trunk/libf/phylmd/surf_land_bucket_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/surf_land_bucket_mod.F90 (revision 3967) +++ LMDZ6/trunk/libf/phylmd/surf_land_bucket_mod.F90 (revision 3974) @@ -24,5 +24,5 @@ USE cpl_mod USE dimphy - USE geometry_mod, ONLY: latitude + USE geometry_mod, ONLY: longitude,latitude USE mod_grid_phy_lmdz USE mod_phys_lmdz_para @@ -103,6 +103,8 @@ ! calculate temperature, heat capacity and conduction flux in soil - IF (soil_model) THEN - CALL soil(dtime, is_ter, knon, snow, tsurf, tsoil, soilcap, soilflux) + IF (soil_model) THEN + CALL soil(dtime, is_ter, knon, snow, tsurf, qsol, & + & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux) + DO i=1, knon cal(i) = RCPD / soilcap(i) Index: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 3967) +++ LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 (revision 3974) @@ -25,4 +25,5 @@ USE dimphy + USE geometry_mod, ONLY : longitude,latitude USE surface_data, ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc USE fonte_neige_mod, ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global @@ -278,5 +279,6 @@ ! use soil model and recalculate properly cal IF (soil_model) THEN - CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux) + CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, & + & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux) cal(1:knon) = RCPD / soilcap(1:knon) radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)