Changeset 2900
- Timestamp:
- Feb 21, 2023, 5:00:24 PM (21 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 9 edited
-
README (modified) (1 diff)
-
libf/phymars/comslope_mod.F90 (modified) (2 diffs)
-
libf/phymars/comsoil_h.F90 (modified) (2 diffs)
-
libf/phymars/dimradmars_mod.F90 (modified) (4 diffs)
-
libf/phymars/phys_state_var_init_mod.F90 (modified) (4 diffs)
-
libf/phymars/physiq_mod.F (modified) (87 diffs)
-
libf/phymars/soil.F (modified) (9 diffs)
-
libf/phymars/surfdat_h.F90 (modified) (3 diffs)
-
libf/phymars/writediagfi.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2897 r2900 3880 3880 Minor edits then (+ svn update with RV had some issues, so there are some "artefact changes" ...) 3881 3881 3882 == 14/02/2023 == RV 3883 First modifications to introduce subslope parametrization in surface processes. 3884 New variables are added in the new module comslope_mod to describe the subslopes: 3885 _ nslope (number of subslope) 3886 _ subslope_dist (distribution of the slopes) 3887 _ def_slope (bound of the slope statistics / repartitions) etc.. 3888 These variables are added to the starfi file. 3889 Other GCM variables, outputs and subroutines are not modified yet. 3890 Only nslope=1 is accepted so far (ie. same as no subslope parametrization) 3891 3882 3892 == 16/02/2023 == RV 3883 3893 MARS PEM : Deep cleaning of variables name and allocate. 3884 3894 All the "dyn to phys" grid change is done in subroutines and not in the main program. 3895 3896 == 21/02/2023 == RV 3897 Following r2896, further implementation of subslope parametrisation. 3898 Carefull ! This is a devolpment revision and it still need improvements and tests. 3899 However, this commit should not change anything for nslope=1. 3900 Only nslope=1 is possible for now! 3901 Utilitaries are not adapted yet. 3902 Dimension of some variables (30 listed below) is changed to add the nslope dimension. 3903 Outputs (diagfi and restartfi) are not changed yet. 3904 nslope is now read and written in the startfi.nc 3905 3906 -
trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90
r2896 r2900 13 13 IMPLICIT NONE 14 14 15 INTEGER, parameter :: nslope = 1! Number of slopes for the statistic (1)15 INTEGER, SAVE :: nslope ! Number of slopes for the statistic (1) 16 16 INTEGER, SAVE :: iflat ! Number of slopes for the statistic (1) 17 17 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope ! bound of the slope statistics / repartitions (°) … … 20 20 REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: subslope_dist ! distribution of the slopes (1) 21 21 INTEGER,SAVE,ALLOCATABLE,DIMENSION(:) :: major_slope ! Index of the subslope that occupies most of themesh (1) 22 !$OMP THREADPRIVATE( def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)22 !$OMP THREADPRIVATE(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope) 23 23 24 24 CONTAINS -
trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
r2887 r2900 17 17 18 18 ! variables (FC: built in firstcall in soil.F) 19 REAL,SAVE,ALLOCATABLE :: tsoil(:,: ) ! sub-surface temperatures (K)20 real,save,allocatable :: mthermdiff(:,: ) ! (FC) mid-layer thermal diffusivity21 real,save,allocatable :: thermdiff(:,: ) ! (FC) inter-layer thermal diffusivity19 REAL,SAVE,ALLOCATABLE :: tsoil(:,:,:) ! sub-surface temperatures (K) 20 real,save,allocatable :: mthermdiff(:,:,:) ! (FC) mid-layer thermal diffusivity 21 real,save,allocatable :: thermdiff(:,:,:) ! (FC) inter-layer thermal diffusivity 22 22 real,save,allocatable :: coefq(:) ! (FC) q_{k+1/2} coefficients 23 real,save,allocatable :: coefd(:,: ) ! (FC) d_k coefficients24 real,save,allocatable :: alph(:,: ) ! (FC) alpha_k coefficients25 real,save,allocatable :: beta(:,: ) ! beta_k coefficients23 real,save,allocatable :: coefd(:,:,:) ! (FC) d_k coefficients 24 real,save,allocatable :: alph(:,:,:) ! (FC) alpha_k coefficients 25 real,save,allocatable :: beta(:,:,:) ! beta_k coefficients 26 26 real,save :: mu 27 27 real,save,allocatable :: flux_geo(:) ! Geothermal Flux (W/m^2) … … 31 31 contains 32 32 33 subroutine ini_comsoil_h(ngrid )33 subroutine ini_comsoil_h(ngrid,nslope) 34 34 35 35 implicit none 36 36 integer,intent(in) :: ngrid ! number of atmospheric columns 37 37 integer,intent(in) :: nslope ! number of sub grid slopes 38 38 allocate(layer(nsoilmx)) !soil layer depths 39 39 allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths 40 40 allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia 41 41 42 allocate(tsoil(ngrid,nsoilmx )) ! soil temperatures42 allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures 43 43 44 allocate(mthermdiff(ngrid,0:nsoilmx-1 ))45 allocate(thermdiff(ngrid,nsoilmx-1 ))44 allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope)) 45 allocate(thermdiff(ngrid,nsoilmx-1,nslope)) 46 46 allocate(coefq(0:nsoilmx-1)) 47 allocate(coefd(ngrid,nsoilmx-1 ))48 allocate(alph(ngrid,nsoilmx-1 ))49 allocate(beta(ngrid,nsoilmx-1 ))47 allocate(coefd(ngrid,nsoilmx-1,nslope)) 48 allocate(alph(ngrid,nsoilmx-1,nslope)) 49 allocate(beta(ngrid,nsoilmx-1,nslope)) 50 50 allocate(flux_geo(ngrid)) 51 51 -
trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90
r2678 r2900 117 117 !! 118 118 REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1) 119 REAL,SAVE,ALLOCATABLE :: fluxrad_sky(: ) ! rad. flux from sky absorbed by surface (W.m-2)120 REAL,SAVE,ALLOCATABLE :: fluxrad(: ) ! Net radiative surface flux (W.m-2)121 REAL,SAVE,ALLOCATABLE :: albedo(:,: ) ! Surface albedo in each solar band119 REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:,:) ! rad. flux from sky absorbed by surface (W.m-2) 120 REAL,SAVE,ALLOCATABLE :: fluxrad(:,:) ! Net radiative surface flux (W.m-2) 121 REAL,SAVE,ALLOCATABLE :: albedo(:,:,:) ! Surface albedo in each solar band 122 122 REAL,SAVE,ALLOCATABLE :: totcloudfrac(:) ! total cloud fraction over the column 123 123 ! aerosol (dust or ice) extinction optical depth at reference wavelength … … 181 181 contains 182 182 183 subroutine ini_dimradmars_mod(ngrid,nlayer )183 subroutine ini_dimradmars_mod(ngrid,nlayer,nslope) 184 184 185 185 implicit none … … 187 187 integer,intent(in) :: ngrid ! number of atmospheric columns 188 188 integer,intent(in) :: nlayer ! number of atmospheric layers 189 189 integer,intent(in) :: nslope ! number of subgrid scale slopes 190 190 nflev=nlayer 191 191 ! ndomainsz=ngrid … … 195 195 ndlo2=ndlon 196 196 197 allocate(albedo(ngrid,2 ))197 allocate(albedo(ngrid,2,nslope)) 198 198 allocate(dtrad(ngrid,nlayer)) 199 allocate(fluxrad_sky(ngrid ))200 allocate(fluxrad(ngrid ))199 allocate(fluxrad_sky(ngrid,nslope)) 200 allocate(fluxrad(ngrid,nslope)) 201 201 allocate(nueffdust(ngrid,nlayer)) 202 202 allocate(totcloudfrac(ngrid)) -
trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
r2628 r2900 61 61 use dust_rad_adjust_mod, only: ini_dust_rad_adjust_mod, & 62 62 end_dust_rad_adjust_mod 63 use comslope_mod, ONLY: nslope 64 use netcdf 65 USE mod_phys_lmdz_para, ONLY: is_master,bcast 66 63 67 IMPLICIT NONE 64 68 … … 71 75 INTEGER,INTENT(in) :: dyn_nqperes 72 76 INTEGER,INTENT(in) :: dyn_nqfils(nq) 77 character(len=10) :: filename ! name of the startfi.nc 78 integer :: ncid, status, nslope_dim_id 79 integer :: nslope_read 80 81 filename="startfi.nc" 82 if(is_master) then 83 status = nf90_open(filename, nf90_nowrite, ncid) 84 if (status /= nf90_noerr) then 85 nslope=1 86 else 87 status = nf90_inq_dimid(ncid, "nslope", nslope_dim_id) 88 if (status /= nf90_noerr) then 89 nslope=1 90 else 91 status = nf90_inquire_dimension(ncid, nslope_dim_id, len = nslope_read) 92 if (status /= nf90_noerr) then 93 call abort_physic("phys_state_var_init","nslope present but not readable",1) 94 else 95 nslope=nslope_read 96 endif 97 endif 98 endif 99 endif 100 call bcast(nslope) 73 101 74 102 ! set dimension and allocate arrays in tracer_mod 75 103 call end_tracer_mod 76 104 call ini_tracer_mod(nq,tname,dyn_nqperes,dyn_nqfils)! MVals: variables isotopes 77 78 105 79 106 ! initialize "print_control" constants/flags ("prt_level","lunout", etc.) … … 101 128 ! allocate "surfdat_h" arrays 102 129 call end_surfdat_h 103 call ini_surfdat_h(ngrid,nq )130 call ini_surfdat_h(ngrid,nq,nslope) 104 131 105 132 ! allocate "comgeomfi_h" arrays … … 109 136 ! allocate "comsoil_h" arrays 110 137 call end_comsoil_h 111 call ini_comsoil_h(ngrid )138 call ini_comsoil_h(ngrid,nslope) 112 139 113 140 ! set some variables in "dimradmars_mod" 114 141 call end_dimradmars_mod 115 call ini_dimradmars_mod(ngrid,nlayer )142 call ini_dimradmars_mod(ngrid,nlayer,nslope) 116 143 117 144 ! allocate arrays in "yomlw_h" -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r2896 r2900 271 271 real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the 272 272 ! size distribution 273 REAL inertiesoil(ngrid,nsoilmx ) ! Time varying subsurface273 REAL inertiesoil(ngrid,nsoilmx,nslope) ! Time varying subsurface 274 274 ! thermal inertia (J.s-1/2.m-2.K-1) 275 275 ! (used only when tifeedback=.true.) … … 301 301 INTEGER l,ig,ierr,igout,iq,tapphys,isoil 302 302 303 REAL fluxsurf_lw(ngrid ) !incident LW (IR) surface flux (W.m-2)304 REAL fluxsurf_dn_sw(ngrid,2 ) ! Incident SW (solar) surface flux (W.m-2)303 REAL fluxsurf_lw(ngrid,nslope) !incident LW (IR) surface flux (W.m-2) 304 REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2) 305 305 REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2) 306 306 REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2) … … 333 333 334 334 c Tendancies due to various processes: 335 REAL dqsurf(ngrid,nq ) ! tendency for tracers on surface (Kg/m2/s)335 REAL dqsurf(ngrid,nq,nslope) ! tendency for tracers on surface (Kg/m2/s) 336 336 REAL zdtlw(ngrid,nlayer) ! (K/s) 337 337 REAL zdtsw(ngrid,nlayer) ! (K/s) … … 340 340 REAL zdtnirco2(ngrid,nlayer) ! (K/s) 341 341 REAL zdtnlte(ngrid,nlayer) ! (K/s) 342 REAL zdtsurf(ngrid ) ! (K/s)342 REAL zdtsurf(ngrid,nslope) ! (K/s) 343 343 REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer) 344 344 REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2) 345 REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid ) ! (K/s)345 REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope) ! (K/s) 346 346 REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2) 347 347 REAL zdhadj(ngrid,nlayer) ! (K/s) 348 348 REAL zdtgw(ngrid,nlayer) ! (K/s) 349 349 REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2) 350 REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid )350 REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope) 351 351 REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) 352 352 353 REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq )353 REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope) 354 354 REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq) 355 355 REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq) … … 357 357 REAL zdqc(ngrid,nlayer,nq) 358 358 REAL zdqcloudco2(ngrid,nlayer,nq) 359 REAL zdqsc(ngrid,nq )359 REAL zdqsc(ngrid,nq,nslope) 360 360 361 361 REAL zdteuv(ngrid,nlayer) ! (K/s) … … 366 366 real*8 PhiEscH,PhiEscH2,PhiEscD 367 367 368 REAL dwatercap(ngrid ), dwatercap_dif(ngrid) ! (kg/m-2)368 REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope) ! (kg/m-2) 369 369 370 370 c Local variable for local intermediate calcul: 371 REAL zflubid(ngrid )371 REAL zflubid(ngrid,nslope) 372 372 REAL zplanck(ngrid),zpopsk(ngrid,nlayer) 373 373 REAL zdum1(ngrid,nlayer) … … 393 393 394 394 REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid) 395 REAL fluxsurf_dn_sw_tot(ngrid ), fluxsurf_up_sw_tot(ngrid)395 REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid) 396 396 character*2 str2 397 397 ! character*5 str5 … … 534 534 535 535 c Sub-grid scale slopes 536 real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K] 537 real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1] 538 real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1] 539 real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2] 540 real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie 541 real :: inertiedat_tmp(ngrid,nsoilmx,nslope) ! temporary inertiedat for soil.F 536 542 integer :: islope 543 real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting 537 544 538 545 logical :: write_restart … … 557 564 558 565 #ifndef MESOSCALE 559 fluxrad(: )=0566 fluxrad(:,:)=0 560 567 wstar(:)=0. 561 568 #endif … … 600 607 PRINT*,'corresponding criterium = ',def_slope_mean(iflat) 601 608 609 DO islope = 1,nslope 610 inertiedat_tmp(:,:,islope) = inertiedat(:,:) 611 ENDDO 612 602 613 #else 603 614 ! MESOSCALE. Supposedly everything is already set in modules. … … 606 617 print*,"check: rad,cpp,g,r,rcp,daysec" 607 618 print*,rad,cpp,g,r,rcp,daysec 608 PRINT*,'check: tsurf ',tsurf(1 ),tsurf(ngrid)609 PRINT*,'check: tsoil ',tsoil(1,1 ),tsoil(ngrid,nsoilmx)619 PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:) 620 PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:) 610 621 PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx) 611 622 PRINT*,'check: midlayer,layer ', mlayer(:),layer(:) 612 623 PRINT*,'check: tracernames ', noms 613 PRINT*,'check: emis ',emis(1 ),emis(ngrid)624 PRINT*,'check: emis ',emis(1,:),emis(ngrid,:) 614 625 PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1) 615 PRINT*,'check: qsurf ',qsurf(1,1 ),qsurf(ngrid,nq)616 PRINT*,'check: co2ice ',qsurf(1,igcm_co2 ),qsurf(ngrid,igcm_co2)626 PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:) 627 PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:) 617 628 !!! 618 629 day_ini = pday … … 645 656 endif 646 657 ! additionnal "academic" initialization of physics 647 write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!" 648 tsurf(:)=pt(:,1) 658 do islope = 1,nslope 659 tsurf(:,islope)=pt(:,1) 660 enddo 649 661 write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" 650 662 do isoil=1,nsoilmx 651 tsoil(1:ngrid,isoil )=tsurf(1:ngrid)663 tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:) 652 664 enddo 653 665 write(*,*) "Physiq: initializing inertiedat !!" … … 683 695 c Thermal inertia feedback: 684 696 IF (tifeedback) THEN 685 CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil) 697 DO islope = 1,nslope 698 CALL soil_tifeedback(ngrid,nsoilmx, 699 s qsurf(:,:,islope),inertiesoil(:,:,islope)) 700 ENDDO 686 701 CALL soil(ngrid,nsoilmx,firstcall,inertiesoil, 687 702 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 688 703 ELSE 689 CALL soil(ngrid,nsoilmx,firstcall,inertiedat ,704 CALL soil(ngrid,nsoilmx,firstcall,inertiedat_tmp, 690 705 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 691 706 ENDIF ! of IF (tifeedback) … … 694 709 & 'PHYSIQ WARNING! Thermal conduction in the soil turned off' 695 710 DO ig=1,ngrid 696 capcal(ig )=1.e5697 fluxgrd(ig )=0.711 capcal(ig,:)=1.e5 712 fluxgrd(ig,:)=0. 698 713 ENDDO 699 714 ENDIF … … 793 808 pdq(:,:,:)=0 794 809 pdpsrf(:)=0 795 zflubid(: )=0796 zdtsurf(: )=0797 dqsurf(:,: )=0810 zflubid(:,:)=0 811 zdtsurf(:,:)=0 812 dqsurf(:,:,:)=0 798 813 dsodust(:,:)=0. 799 814 dsords(:,:)=0. 800 815 dsotop(:,:)=0. 801 dwatercap(: )=0816 dwatercap(:,:)=0 802 817 803 818 ! Dust scenario conversion coefficient from IRabs to VISext … … 899 914 & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) 900 915 end do 901 co2totA = co2totA + qsurf(ig,igcm_co2) 916 do islope = 1,nslope 917 co2totA = co2totA + qsurf(ig,igcm_co2,islope)* 918 & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 919 enddo 902 920 end do 903 921 else … … 909 927 & +pdq(ig,l,igcm_co2)*ptimestep) 910 928 end do 911 co2totA = co2totA + qsurf(ig,igcm_co2) 929 do islope = 1,nslope 930 co2totA = co2totA + qsurf(ig,igcm_co2,islope)* 931 & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 932 enddo 912 933 end do 913 934 endif ! of if (igcm_co2_ice.ne.0) … … 1007 1028 ! callradite for the part with clouds 1008 1029 clearsky=.false. ! part with clouds for both cases CLFvarying true/false 1009 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, 1010 & emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout, 1011 & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_dn_sw,fluxsurf_up_sw, 1030 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, 1031 & albedo(:,:,1), 1032 & emis(:,1),mu0,zplev,zplay,pt,tsurf(:,1),fract,dist_sol,igout, 1033 & zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat), 1034 & fluxsurf_up_sw, 1012 1035 & fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw, 1013 1036 & tau_pref_scenario,tau_pref_gcm, 1014 1037 & tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, 1015 1038 & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2, 1016 & qsurf(:,igcm_co2 ),rstormdust,rtopdust,totstormfract,1039 & qsurf(:,igcm_co2,1),rstormdust,rtopdust,totstormfract, 1017 1040 & clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac) 1041 1042 DO islope=1,nslope 1043 fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat) 1044 fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat) 1045 ENDDO 1018 1046 1019 1047 ! case of sub-grid water ice clouds: callradite for the clear case … … 1025 1053 clearsky=.true. 1026 1054 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, 1027 & albedo,emis,mu0,zplev,zplay,pt,tsurf,fract, 1055 & albedo(:,:,1),emis(:,1),mu0,zplev,zplay,pt, 1056 & tsurf(:,1),fract, 1028 1057 & dist_sol,igout,zdtlwclf,zdtswclf, 1029 1058 & fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf, … … 1032 1061 & dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, 1033 1062 & taucloudtesclf,rdust, 1034 & rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2 ),1063 & rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2,1), 1035 1064 & rstormdust,rtopdust,totstormfract, 1036 1065 & clearatm,dsords,dsotop, … … 1042 1071 tf_clf=totcloudfrac(ig) 1043 1072 ntf_clf=1.-tf_clf 1044 fluxsurf_lw(ig) = ntf_clf*fluxsurf_lwclf(ig) 1045 & + tf_clf*fluxsurf_lw(ig) 1046 fluxsurf_dn_sw(ig,1:2) = 1073 DO islope=1,nslope 1074 fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) 1075 & + tf_clf*fluxsurf_lw(ig,islope) 1076 fluxsurf_dn_sw(ig,1:2,islope) = 1047 1077 & ntf_clf*fluxsurf_dn_swclf(ig,1:2) 1048 & + tf_clf*fluxsurf_dn_sw(ig,1:2) 1078 & + tf_clf*fluxsurf_dn_sw(ig,1:2,islope) 1079 ENDDO 1049 1080 fluxsurf_up_sw(ig,1:2) = 1050 1081 & ntf_clf*fluxsurf_up_swclf(ig,1:2) … … 1099 1130 c --------------------------------------------------------- 1100 1131 IF(callslope) THEN 1132 ! assume that in this case, nslope = 1 1133 if(nslope.ne.1) then 1134 call abort_physic( 1135 & "physiq","callslope=true but nslope.ne.1",1) 1136 endif 1101 1137 print *, 'Slope scheme is on and computing...' 1102 1138 DO ig=1,ngrid 1103 1139 sl_the = theta_sl(ig) 1104 1140 IF (sl_the .ne. 0.) THEN 1105 ztim1=fluxsurf_dn_sw(ig,1 )+fluxsurf_dn_sw(ig,2)1141 ztim1=fluxsurf_dn_sw(ig,1,1)+fluxsurf_dn_sw(ig,2,1) 1106 1142 DO l=1,2 1107 1143 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. … … 1109 1145 sl_lat = 180.*latitude(ig)/pi 1110 1146 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) 1111 sl_alb = albedo(ig,l )1147 sl_alb = albedo(ig,l,1) 1112 1148 sl_psi = psi_sl(ig) 1113 sl_fl0 = fluxsurf_dn_sw(ig,l )1149 sl_fl0 = fluxsurf_dn_sw(ig,l,1) 1114 1150 sl_di0 = 0. 1115 1151 if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then … … 1117 1153 sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol 1118 1154 sl_di0 = sl_di0/ztim1 1119 sl_di0 = fluxsurf_dn_sw(ig,l )*sl_di01155 sl_di0 = fluxsurf_dn_sw(ig,l,1)*sl_di0 1120 1156 endif 1121 1157 ! you never know (roundup concern...) … … 1126 1162 & sl_di0, sl_fl0, sl_flu ) 1127 1163 !!!!!!!!!!!!!!!!!!!!!!!!!! 1128 fluxsurf_dn_sw(ig,l ) = sl_flu1164 fluxsurf_dn_sw(ig,l,1) = sl_flu 1129 1165 ENDDO 1130 1166 !!! compute correction on IR flux as well 1131 1167 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 1132 fluxsurf_lw(ig )= fluxsurf_lw(ig)*sky1168 fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky 1133 1169 ENDIF 1134 1170 ENDDO 1135 ENDIF 1171 c ELSE ! not calling subslope, nslope might be > 1 1172 c DO islope = 1,nslope 1173 c IF(def_slope_mean(islope).ge.0.) THEN 1174 c sl_psi = 0. !Northward slope 1175 c ELSE 1176 c sl_psi = 180. !Southward slope 1177 c ENDIF 1178 c sl_the=abs(def_slope_mean(islope)) 1179 c IF (sl_the .gt. 1e-6) THEN 1180 c DO ig=1,ngrid 1181 c ztim1=fluxsurf_dn_sw(ig,1,islope) 1182 c s +fluxsurf_dn_sw(ig,2,islope) 1183 c DO l=1,2 1184 c sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. 1185 c sl_ra = pi*(1.0-sl_lct/12.) 1186 c sl_lat = 180.*latitude(ig)/pi 1187 c sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) 1188 c sl_alb = albedo(ig,l,islope) 1189 c sl_psi = psi_sl(ig) 1190 c sl_fl0 = fluxsurf_dn_sw(ig,l,islope) 1191 c sl_di0 = 0. 1192 c if (mu0(ig) .gt. 0.) then 1193 c sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) 1194 c sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol 1195 c sl_di0 = sl_di0/ztim1 1196 c sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0 1197 c endif 1198 c ! you never know (roundup concern...) 1199 c if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 1200 c !!!!!!!!!!!!!!!!!!!!!!!!!! 1201 c CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 1202 c & sl_tau, sl_alb, sl_the, sl_psi, 1203 c & sl_di0, sl_fl0, sl_flu ) 1204 c !!!!!!!!!!!!!!!!!!!!!!!!!! 1205 c fluxsurf_dn_sw(ig,l,islope) = sl_flu 1206 c ENDDO 1207 c !!! compute correction on IR flux as well 1208 c 1209 c fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope) 1210 c & *sky_slope(islope) 1211 c ENDDO 1212 c ENDIF 1213 c ENDDO ! islope = 1,nslope 1214 ENDIF ! callslope 1136 1215 1137 1216 c CO2 near infrared absorption … … 1146 1225 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1147 1226 DO ig=1,ngrid 1148 fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) 1149 $ +fluxsurf_dn_sw(ig,1)*(1.-albedo(ig,1)) 1150 $ +fluxsurf_dn_sw(ig,2)*(1.-albedo(ig,2)) 1227 DO islope = 1,nslope 1228 fluxrad_sky(ig,islope) = 1229 $ emis(ig,islope)*fluxsurf_lw(ig,islope) 1230 $ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope)) 1231 $ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope)) 1232 ENDDO 1151 1233 ENDDO 1152 1234 … … 1175 1257 c 1176 1258 DO ig=1,ngrid 1177 zplanck(ig)=tsurf(ig)*tsurf(ig) 1178 zplanck(ig)=emis(ig)* 1259 DO islope = 1,nslope 1260 zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope) 1261 zplanck(ig)=emis(ig,islope)* 1179 1262 $ stephan*zplanck(ig)*zplanck(ig) 1180 fluxrad(ig )=fluxrad_sky(ig)-zplanck(ig)1263 fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig) 1181 1264 IF(callslope) THEN 1182 1265 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 1183 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) 1266 fluxrad(ig,nslope)=fluxrad(ig,nslope)+ 1267 $ (1.-sky)*zplanck(ig) 1268 ELSE 1269 fluxrad(ig,islope)=fluxrad(ig,islope) + 1270 $ (1.-sky_slope(iflat))*emis(ig,iflat)* 1271 $ stephan*tsurf(ig,iflat)**4 1184 1272 ENDIF 1185 1273 ENDDO 1274 ENDDO 1186 1275 1187 1276 DO l=1,nlayer … … 1224 1313 c for radiative transfer 1225 1314 & clearatm,icount,zday,zls, 1226 & tsurf,qsurf(:,igcm_co2 ),igout,1315 & tsurf,qsurf(:,igcm_co2,iflat),igout, 1227 1316 & totstormfract,tauscaling, 1228 1317 & dust_rad_adjust,IRtoVIScoef, … … 1287 1376 & pq,pdq,pt,pdt,zplev,zplay,zzlev, 1288 1377 & zzlay,zdtsw,zdtlw, 1289 & icount,zday,zls,tsurf,qsurf(:,igcm_co2), 1378 & icount,zday,zls,tsurf(:,iflat), 1379 & qsurf(:,igcm_co2,iflat), 1290 1380 & igout,aerosol,tauscaling,dust_rad_adjust, 1291 1381 & IRtoVIScoef,totstormfract,clearatm, … … 1349 1439 IF (calldifv) THEN 1350 1440 DO ig=1,ngrid 1351 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) 1441 DO islope = 1,nslope 1442 zflubid(ig,islope)=fluxrad(ig,islope) 1443 & +fluxgrd(ig,islope) 1444 ENDDO 1352 1445 ENDDO 1353 1354 1446 zdum1(:,:)=0 1355 1447 zdum2(:,:)=0 … … 1371 1463 1372 1464 DO ig=1, ngrid 1373 IF (zh(ig,1) .lt. tsurf(ig )) THEN1465 IF (zh(ig,1) .lt. tsurf(ig,1)) THEN 1374 1466 wstar(ig)=1. 1375 1467 hfmax_th(ig)=0.2 … … 1390 1482 1391 1483 c Calling vdif (Martian version WITH CO2 condensation) 1392 dwatercap_dif(: ) = 0.1484 dwatercap_dif(:,:) = 0. 1393 1485 zcdh(:) = 0. 1394 1486 zcdv(:) = 0. … … 1405 1497 1406 1498 DO ig=1,ngrid 1407 zdtsurf(ig )=zdtsurf(ig)+zdtsdif(ig)1408 dwatercap(ig )=dwatercap(ig)+dwatercap_dif(ig)1499 zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:) 1500 dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:) 1409 1501 ENDDO 1410 1502 c call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif, 1503 c & albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp) 1411 1504 IF (.not.turb_resolved) THEN 1412 1505 DO l=1,nlayer … … 1429 1522 DO iq=1, nq 1430 1523 DO ig=1,ngrid 1431 dqsurf(ig,iq )=dqsurf(ig,iq) + zdqsdif(ig,iq)1524 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) 1432 1525 ENDDO 1433 1526 ENDDO … … 1446 1539 zdqdif(1:ngrid,1,1:nq)=0. 1447 1540 zdqdif(1:ngrid,1,igcm_dust_number) = 1448 . -zdqsdif(1:ngrid,igcm_dust_number )1541 . -zdqsdif(1:ngrid,igcm_dust_number,iflat) 1449 1542 zdqdif(1:ngrid,1,igcm_dust_mass) = 1450 . -zdqsdif(1:ngrid,igcm_dust_mass )1543 . -zdqsdif(1:ngrid,igcm_dust_mass,iflat) 1451 1544 zdqdif(1:ngrid,2:nlayer,1:nq) = 0. 1452 1545 DO iq=1, nq 1453 1546 IF ((iq .ne. igcm_dust_mass) 1454 1547 & .and. (iq .ne. igcm_dust_number)) THEN 1455 zdqsdif(:,iq )=0.1548 zdqsdif(:,iq,:)=0. 1456 1549 ENDIF 1457 1550 ENDDO 1458 1551 ELSE 1459 1552 zdqdif(1:ngrid,1:nlayer,1:nq) = 0. 1460 zdqsdif(1:ngrid,1:nq ) = 0.1553 zdqsdif(1:ngrid,1:nq,1:nslope) = 0. 1461 1554 ENDIF 1462 1555 ENDIF 1463 1556 ELSE 1464 1557 DO ig=1,ngrid 1465 zdtsurf(ig)=zdtsurf(ig)+ 1466 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 1558 DO islope=1,nslope 1559 zdtsurf(ig,islope)=zdtsurf(ig,islope)+ 1560 s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope) 1561 ENDDO 1467 1562 ENDDO 1468 1563 IF (turb_resolved) THEN … … 1707 1802 DO iq=1, nq 1708 1803 DO ig=1,ngrid 1709 dqsurf(ig,iq)=dqsurf(ig,iq)+zdqssed_ccn(ig,iq) 1804 DO islope = 1,nslope 1805 dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+ 1806 & zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.) 1807 ENDDO !(islope) 1710 1808 ENDDO ! (ig) 1711 ENDDO ! (iq) 1809 ENDDO ! (iq)q) 1712 1810 c Temperature variation due to latent heat release 1713 1811 pdt(1:ngrid,1:nlayer) = … … 1884 1982 do iq=1,nq 1885 1983 DO ig=1,ngrid 1886 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq) 1984 DO islope = 1,nslope 1985 dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + 1986 & zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.) 1987 ENDDO 1887 1988 ENDDO 1888 1989 enddo … … 1933 2034 DO iq=1, nq 1934 2035 DO ig=1,ngrid 1935 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq) 2036 DO islope = 1,nslope 2037 dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + 2038 & zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.) 2039 ENDDO 1936 2040 ENDDO 1937 2041 ENDDO … … 1950 2054 DO iq=1, nq 1951 2055 DO ig=1,ngrid 1952 dqsurf(ig,iq )=dqsurf(ig,iq) + zdqsdif(ig,iq)2056 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) 1953 2057 ENDDO 1954 2058 ENDDO … … 1974 2078 $ surfdust, surfice) 1975 2079 ! call photochemistry 2080 DO ig = 1,ngrid 2081 qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig)) 2082 ENDDO 1976 2083 call calchim(ngrid,nlayer,nq, 1977 2084 & ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0, 1978 2085 $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, 1979 $ zdqcloud,zdqscloud,tau(:,1),qsurf(:,igcm_co2), 2086 $ zdqcloud,zdqscloud,tau(:,1), 2087 $ qsurf_tmp(:,igcm_co2), 1980 2088 $ pu,pdu,pv,pdv,surfdust,surfice) 1981 2089 … … 2006 2114 ! tracers is zero anyways 2007 2115 DO ig=1,ngrid 2008 dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) 2116 DO islope = 1,nslope 2117 dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) + 2118 & zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.) 2119 ENDDO 2009 2120 ENDDO 2010 2121 ENDDO ! of DO iq=1,nq … … 2013 2124 if (igcm_h2o2.ne.0) then 2014 2125 DO ig=1,ngrid 2015 dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) 2016 & +zdqscloud(ig,igcm_h2o2) 2126 DO islope = 1,nslope 2127 dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+ 2128 & zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.) 2129 ENDDO 2017 2130 ENDDO 2018 2131 endif … … 2029 2142 if (callthermos) then 2030 2143 call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol, 2031 $ mu0,ptimestep,ptime,zday,tsurf ,zzlev,zzlay,2144 $ mu0,ptimestep,ptime,zday,tsurf(:,iflat),zzlev,zzlay, 2032 2145 & pt,pq,pu,pv,pdt,pdq, 2033 2146 $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff, … … 2056 2169 CALL geticecover(ngrid, 180.*zls/pi, 2057 2170 . 180.*longitude/pi, 180.*latitude/pi, 2058 . qsurf (:,igcm_co2) )2059 qsurf (:,igcm_co2) = qsurf(:,igcm_co2) * 10000.2171 . qsurf_tmp(:,igcm_co2) ) 2172 qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000. 2060 2173 ENDIF 2061 2174 … … 2063 2176 IF (callcond) THEN 2064 2177 zdtc(:,:) = 0. 2065 zdtsurfc(: ) = 0.2178 zdtsurfc(:,:) = 0. 2066 2179 zduc(:,:) = 0. 2067 2180 zdvc(:,:) = 0. 2068 2181 zdqc(:,:,:) = 0. 2069 zdqsc(:,: ) = 0.2182 zdqsc(:,:,:) = 0. 2070 2183 CALL co2condens(ngrid,nlayer,nq,ptimestep, 2071 2184 $ capcal,zplay,zplev,tsurf,pt, 2072 2185 $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, 2073 $ qsurf(:,igcm_co2 ),albedo,emis,rdust,2186 $ qsurf(:,igcm_co2,:),albedo,emis,rdust, 2074 2187 $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, 2075 2188 $ fluxsurf_dn_sw,zls, … … 2079 2192 DO iq=1, nq 2080 2193 DO ig=1,ngrid 2081 dqsurf(ig,iq )=dqsurf(ig,iq)+zdqsc(ig,iq)2194 dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:) 2082 2195 ENDDO ! (ig) 2083 2196 ENDDO ! (iq) … … 2090 2203 ENDDO 2091 2204 DO ig=1,ngrid 2092 zdtsurf(ig ) = zdtsurf(ig) + zdtsurfc(ig)2205 zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:) 2093 2206 ENDDO 2094 2207 … … 2130 2243 DO iq=1, nq 2131 2244 DO ig=1,ngrid 2132 2133 qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq) 2134 2245 DO islope = 1,nslope 2246 qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+ 2247 & ptimestep*dqsurf(ig,iq,islope) 2248 ENDDO 2135 2249 ENDDO ! (ig) 2136 2250 ENDDO ! (iq) … … 2144 2258 2145 2259 DO ig=1,ngrid 2146 tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 2260 DO islope = 1,nslope 2261 tsurf(ig,islope)=tsurf(ig,islope)+ 2262 & ptimestep*zdtsurf(ig,islope) 2263 ENDDO 2147 2264 ENDDO 2148 2265 … … 2169 2286 c ------------------------------------------------------------- 2170 2287 do ig=1,ngrid 2171 if ((qsurf(ig,igcm_co2).eq.0).and. 2172 & (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then 2288 do islope = 1,nslope 2289 if ((qsurf(ig,igcm_co2,islope).eq.0).and. 2290 & (qsurf(ig,igcm_h2o_ice,islope) 2291 & .gt.frost_albedo_threshold)) then 2173 2292 if ((watercaptag(ig)).and.(cst_cap_albedo)) then 2174 albedo(ig,1 ) = albedo_h2o_cap2175 albedo(ig,2 ) = albedo_h2o_cap2293 albedo(ig,1,islope) = albedo_h2o_cap 2294 albedo(ig,2,islope) = albedo_h2o_cap 2176 2295 else 2177 albedo(ig,1 ) = albedo_h2o_frost2178 albedo(ig,2 ) = albedo_h2o_frost2296 albedo(ig,1,islope) = albedo_h2o_frost 2297 albedo(ig,2,islope) = albedo_h2o_frost 2179 2298 endif !((watercaptag(ig)).and.(cst_cap_albedo)) then 2180 2299 c write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice) … … 2182 2301 c & ,latitude(ig)*180./pi, longitude(ig)*180./pi 2183 2302 endif 2303 enddo ! islope 2184 2304 enddo ! of do ig=1,ngrid 2185 2305 ENDIF ! of IF (water) … … 2191 2311 c Thermal inertia feedback 2192 2312 IF (tifeedback) THEN 2193 CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil) 2313 DO islope = 1,nslope 2314 CALL soil_tifeedback(ngrid,nsoilmx, 2315 s qsurf(:,:,islope),inertiesoil(:,:,islope)) 2316 ENDDO 2194 2317 CALL soil(ngrid,nsoilmx,.false.,inertiesoil, 2195 2318 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 2196 2319 ELSE 2197 CALL soil(ngrid,nsoilmx,.false.,inertiedat ,2320 CALL soil(ngrid,nsoilmx,.false.,inertiedat_tmp, 2198 2321 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 2199 2322 ENDIF … … 2242 2365 2243 2366 DO ig=1,ngrid 2244 watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig) 2367 DO islope = 1,nslope 2368 watercap(ig,islope)=watercap(ig,islope)+ 2369 s ptimestep*dwatercap(ig,islope) 2370 ENDDO 2245 2371 ENDDO 2246 2372 … … 2248 2374 2249 2375 DO ig=1,ngrid 2250 if (watercaptag(ig).and. 2251 & (qsurf(ig,igcm_h2o_ice).gt.frost_metam_threshold)) then 2252 2253 watercap(ig)=watercap(ig)+qsurf(ig,igcm_h2o_ice) 2254 & - frost_metam_threshold 2255 qsurf(ig,igcm_h2o_ice) = frost_metam_threshold 2376 DO islope = 1,nslope 2377 if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope) 2378 & .gt.frost_metam_threshold)) then 2379 2380 watercap(ig,islope)=watercap(ig,islope) 2381 & +qsurf(ig,igcm_h2o_ice,islope) 2382 & - frost_metam_threshold 2383 qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold 2256 2384 endif ! (watercaptag(ig).and. 2385 ENDDO 2257 2386 ENDDO 2258 2387 … … 2263 2392 c 13. Write output files 2264 2393 c ---------------------- 2394 2395 call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf, 2396 & albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg) 2265 2397 2266 2398 c ------------------------------- … … 2312 2444 fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) + 2313 2445 & fluxtop_up_sw(1:ngrid,2) 2314 fluxsurf_dn_sw_tot(1:ngrid)=fluxsurf_dn_sw(1:ngrid,1) + 2315 & fluxsurf_dn_sw(1:ngrid,2) 2446 fluxsurf_dn_sw_tot(1:ngrid,1:nslope)= 2447 & fluxsurf_dn_sw(1:ngrid,1,1:nslope) + 2448 & fluxsurf_dn_sw(1:ngrid,2,1:nslope) 2316 2449 fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) + 2317 2450 & fluxsurf_up_sw(1:ngrid,2) … … 2344 2477 WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' 2345 2478 WRITE(*,'(2f10.5,2f15.5)') 2346 s tsurf(igout ),zdtsurf(igout)*ptimestep,2479 s tsurf(igout,:),zdtsurf(igout,:)*ptimestep, 2347 2480 s zplev(igout,1),pdpsrf(igout)*ptimestep 2348 2481 WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' … … 2374 2507 2375 2508 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 2376 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf ,zh,z_out,n_out,2377 & T_out,u_out,ustar,tstar,L_mo,vhf,vvv)2509 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat), 2510 & zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv) 2378 2511 ! pourquoi ustar recalcule ici? fait dans vdifc. 2379 2512 … … 2472 2605 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, 2473 2606 . ptimestep,ztime_fin, 2474 . tsurf,tsoil,albedo,emis, 2475 . q2,qsurf,tauscaling,totcloudfrac,wstar, 2476 . watercap) 2607 . tsurf(:,iflat),tsoil(:,:,iflat),albedo(:,:,iflat), 2608 . emis(:,iflat), 2609 . q2,qsurf(:,:,iflat),tauscaling,totcloudfrac,wstar, 2610 . watercap(:,iflat)) 2477 2611 2478 2612 ENDIF ! of IF (write_restart) … … 2677 2811 2678 2812 call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) 2679 call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) 2813 call wstats(ngrid,"tsurf","Surface temperature","K",2 2814 & ,tsurf(:,iflat)) 2680 2815 call wstats(ngrid,"co2ice","CO2 ice cover", 2681 & "kg.m-2",2,qsurf(:,igcm_co2 ))2816 & "kg.m-2",2,qsurf(:,igcm_co2,iflat)) 2682 2817 call wstats(ngrid,"watercap","H2O ice cover", 2683 & "kg.m-2",2,watercap )2818 & "kg.m-2",2,watercap(:,iflat)) 2684 2819 call wstats(ngrid,"tau_pref_scenario", 2685 2820 & "prescribed visible dod at 610 Pa","NU", … … 2690 2825 call wstats(ngrid,"fluxsurf_lw", 2691 2826 & "Thermal IR radiative flux to surface","W.m-2",2, 2692 & fluxsurf_lw )2827 & fluxsurf_lw(:,iflat)) 2693 2828 call wstats(ngrid,"fluxsurf_dn_sw", 2694 2829 & "Incoming Solar radiative flux to surface","W.m-2",2, 2695 & fluxsurf_dn_sw_tot )2830 & fluxsurf_dn_sw_tot(:,iflat)) 2696 2831 call wstats(ngrid,"fluxsurf_up_sw", 2697 2832 & "Reflected Solar radiative flux from surface","W.m-2",2, … … 2718 2853 & "m2.s-2",3,q2) 2719 2854 call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 2720 & emis )2855 & emis(:,iflat)) 2721 2856 c call wstats(ngrid,"ssurf","Surface stress","N.m-2", 2722 2857 c & 2,zstress) … … 2764 2899 call wstats(ngrid,"h2o_ice_s", 2765 2900 & "surface h2o_ice","kg/m2", 2766 & 2,qsurf(1,igcm_h2o_ice ))2901 & 2,qsurf(1,igcm_h2o_ice,iflat)) 2767 2902 call wstats(ngrid,'albedo', 2768 2903 & 'albedo', 2769 & '',2,albedo(1,1 ))2904 & '',2,albedo(1,1,iflat)) 2770 2905 call wstats(ngrid,"mtot", 2771 2906 & "total mass of water vapor","kg/m2", … … 2852 2987 & 'kg/kg',3,zq(1,1,iq)) 2853 2988 call wstats(ngrid,'qsurf'//str2,'qsurf', 2854 & 'kg.m-2',2,qsurf(1,iq ))2989 & 'kg.m-2',2,qsurf(1,iq,iflat)) 2855 2990 end do 2856 2991 endif ! (doubleq) … … 3019 3154 c for any variables !! 3020 3155 call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 3021 & emis )3156 & emis(:,iflat)) 3022 3157 c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) 3023 3158 c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) … … 3029 3164 & "m2s-2",2,phisfi) 3030 3165 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 3031 & tsurf )3166 & tsurf(:,iflat)) 3032 3167 call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 3033 3168 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness" 3034 & ,"kg.m-2",2,qsurf(:,igcm_co2 ))3169 & ,"kg.m-2",2,qsurf(:,igcm_co2,iflat)) 3035 3170 call WRITEDIAGFI(ngrid,"watercap","Water ice thickness" 3036 & ,"kg.m-2",2,watercap)3171 & ,"kg.m-2",2,watercap(:,iflat)) 3037 3172 3038 3173 call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1", … … 3041 3176 & "K",2,zt(1,7)) 3042 3177 call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, 3043 & fluxsurf_lw )3178 & fluxsurf_lw(:,iflat)) 3044 3179 call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw","fluxsurf_dn_sw", 3045 & "W.m-2",2,fluxsurf_dn_sw_tot )3180 & "W.m-2",2,fluxsurf_dn_sw_tot(:,iflat)) 3046 3181 call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, 3047 3182 & fluxtop_lw) … … 3079 3214 !!! this is to ensure correct initialisation of mesoscale model 3080 3215 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 3081 & tsurf )3216 & tsurf(:,iflat)) 3082 3217 call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 3083 3218 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, 3084 & qsurf(:,igcm_co2 ))3219 & qsurf(:,igcm_co2,iflat)) 3085 3220 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 3086 3221 call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) … … 3089 3224 & emis) 3090 3225 call WRITEDIAGFI(ngrid,"tsoil","Soil temperature", 3091 & "K",3,tsoil )3226 & "K",3,tsoil(:,:,iflat)) 3092 3227 call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia", 3093 3228 & "K",3,inertiedat) … … 3170 3305 call WRITEDIAGFI(ngrid,'qsurf02','surface tracer', 3171 3306 & 'kg.m-2',2, 3172 & qsurf(1:ngrid,igcm_h2o_ice ))3307 & qsurf(1:ngrid,igcm_h2o_ice,iflat)) 3173 3308 #endif 3174 3309 call WRITEDIAGFI(ngrid,'mtot', … … 3265 3400 call WRITEDIAGFI(ngrid,'h2o_ice_s', 3266 3401 & 'surface h2o_ice', 3267 & 'kg.m-2',2,qsurf(1,igcm_h2o_ice ))3402 & 'kg.m-2',2,qsurf(1,igcm_h2o_ice,iflat)) 3268 3403 if (hdo) then 3269 3404 call WRITEDIAGFI(ngrid,'hdo_ice_s', 3270 3405 & 'surface hdo_ice', 3271 & 'kg.m-2',2,qsurf(1,igcm_hdo_ice ))3406 & 'kg.m-2',2,qsurf(1,igcm_hdo_ice,iflat)) 3272 3407 3273 3408 do ig=1,ngrid 3274 if (qsurf (ig,igcm_h2o_ice).gt.qperemin) then3275 DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/3276 & qsurf(ig,igcm_h2o_ice) )/155.76e-63409 if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then 3410 DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/ 3411 & qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6 3277 3412 else 3278 3413 DoH_surf(ig) = 0. … … 3287 3422 CALL WRITEDIAGFI(ngrid,'albedo', 3288 3423 & 'albedo', 3289 & '',2,albedo(1,1 ))3424 & '',2,albedo(1,1,iflat)) 3290 3425 if (tifeedback) then 3291 3426 call WRITEDIAGSOIL(ngrid,"soiltemp", … … 3294 3429 call WRITEDIAGSOIL(ngrid,'soilti', 3295 3430 & 'Soil Thermal Inertia', 3296 & 'J.s-1/2.m-2.K-1',3,inertiesoil )3431 & 'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,iflat)) 3297 3432 endif 3298 3433 !A. Pottier … … 3392 3527 & 'kg/kg',3,zq(1,1,iq)) 3393 3528 call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf', 3394 & 'kg.m-2',2,qsurf(1,iq ))3529 & 'kg.m-2',2,qsurf(1,iq,iflat)) 3395 3530 end do 3396 3531 endif ! (doubleq) … … 3419 3554 call WRITEDIAGFI(ngrid,'qsurf', 3420 3555 & 'stormdust injection', 3421 & 'kg.m-2',2,qsurf(:,igcm_stormdust_mass ))3556 & 'kg.m-2',2,qsurf(:,igcm_stormdust_mass,iflat)) 3422 3557 call WRITEDIAGFI(ngrid,'pdqsurf', 3423 3558 & 'tendancy stormdust mass at surface', 3424 & 'kg.m-2',2,dqsurf(:,igcm_stormdust_mass ))3559 & 'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,iflat)) 3425 3560 call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust', 3426 3561 & 'm/s',3,wspeed(:,1:nlayer)) … … 3484 3619 & 'part/kg',3,nccn) 3485 3620 call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr', 3486 & 'kg.m-2',2,qsurf(1,igcm_ccn_mass ))3621 & 'kg.m-2',2,qsurf(1,igcm_ccn_mass,iflat)) 3487 3622 call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number', 3488 & 'kg.m-2',2,qsurf(1,igcm_ccn_number ))3623 & 'kg.m-2',2,qsurf(1,igcm_ccn_number,iflat)) 3489 3624 endif ! (scavenging) 3490 3625 … … 3636 3771 ! Write soil temperature 3637 3772 call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 3638 & 3,tsoil )3773 & 3,tsoil(:,:,iflat)) 3639 3774 ! Write surface temperature 3640 3775 ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", … … 3666 3801 c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') 3667 3802 call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 3668 & 3,tsoil )3803 & 3,tsoil(:,:,iflat)) 3669 3804 3670 3805 ! THERMALS STUFF 1D … … 3690 3825 call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) 3691 3826 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, 3692 & tsurf )3827 & tsurf(:,iflat)) 3693 3828 call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu) 3694 3829 call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv) … … 3704 3839 & 'w.m-2',1,zdtlw) 3705 3840 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness" 3706 & ,"kg.m-2",0,qsurf(:,igcm_co2))3841 & ,"kg.m-2",0,qsurf(:,igcm_co2,iflat)) 3707 3842 3708 3843 if (igcm_co2.ne.0) then … … 3742 3877 & ,1,aerosol(:,:,iaer_stormdust_doubleq)) 3743 3878 call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion', 3744 & 'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))3879 & 'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass,iflat)) 3745 3880 call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion', 3746 & 'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))3881 & 'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass,iflat)) 3747 3882 call WRITEDIAGFI(ngrid,'mstormdtot', 3748 3883 & 'total mass of stormdust only', … … 3767 3902 call WRITEDIAGFI(ngrid,'rdsqsurf', 3768 3903 & 'stormdust at surface', 3769 & 'kg.m-2',0,qsurf(:,igcm_stormdust_mass ))3904 & 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,iflat)) 3770 3905 call WRITEDIAGFI(ngrid,'qsurf', 3771 3906 & 'dust mass at surface', 3772 & 'kg.m-2',0,qsurf(:,igcm_dust_mass ))3907 & 'kg.m-2',0,qsurf(:,igcm_dust_mass,iflat)) 3773 3908 call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust', 3774 3909 & 'm/s',1,wspeed) … … 3824 3959 mtot = 0 3825 3960 icetot = 0 3826 h2otot = qsurf(1,igcm_h2o_ice )3961 h2otot = qsurf(1,igcm_h2o_ice,iflat) 3827 3962 if (hdo) THEN 3828 3963 mtotD = 0 3829 3964 icetotD = 0 3830 hdotot = qsurf(1,igcm_hdo_ice )3965 hdotot = qsurf(1,igcm_hdo_ice,iflat) 3831 3966 ENDIF !hdo 3832 3967 … … 3860 3995 CALL WRITEDIAGFI(ngrid,'albedo', 3861 3996 & 'albedo', 3862 & '',2,albedo(1,1 ))3997 & '',2,albedo(1,1,iflat)) 3863 3998 3864 3999 IF (hdo) THEN … … 3970 4105 do iq=1,nq 3971 4106 call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s', 3972 & trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))4107 & trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq,iflat)) 3973 4108 end do 3974 4109 … … 4030 4165 & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) 4031 4166 enddo 4032 co2totB = co2totB + qsurf(ig,igcm_co2 )4167 co2totB = co2totB + qsurf(ig,igcm_co2,iflat) 4033 4168 enddo 4034 4169 else … … 4039 4174 & (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep) 4040 4175 enddo 4041 co2totB = co2totB + qsurf(ig,igcm_co2 )4176 co2totB = co2totB + qsurf(ig,igcm_co2,iflat) 4042 4177 enddo 4043 4178 endif ! of if (igcm_co2_ice.ne.0) -
trunk/LMDZ.MARS/libf/phymars/soil.F
r2887 r2900 8 8 & coefd, alph, beta, mu,flux_geo 9 9 use surfdat_h, only: watercaptag, inert_h2o_ice 10 10 use comslope_mod, ONLY: nslope 11 11 implicit none 12 12 … … 29 29 integer nsoil ! number of soil layers 30 30 logical firstcall ! identifier for initialization call 31 real therm_i(ngrid,nsoil ) ! thermal inertia31 real therm_i(ngrid,nsoil,nslope) ! thermal inertia 32 32 real timestep ! time step 33 real tsurf(ngrid ) ! surface temperature33 real tsurf(ngrid,nslope) ! surface temperature 34 34 ! outputs: 35 real tsoil(ngrid,nsoil ) ! soil (mid-layer) temperature36 real capcal(ngrid ) ! surface specific heat37 real fluxgrd(ngrid ) ! surface diffusive heat flux35 real tsoil(ngrid,nsoil,nslope) ! soil (mid-layer) temperature 36 real capcal(ngrid,nslope) ! surface specific heat 37 real fluxgrd(ngrid,nslope) ! surface diffusive heat flux 38 38 39 39 ! local variables: 40 integer ig,ik 40 integer ig,ik,islope 41 41 42 42 ! 0. Initialisations and preprocessing step … … 46 46 ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities 47 47 do ig=1,ngrid 48 do islope = 1,nslope 48 49 if (watercaptag(ig)) then 49 50 do ik=0,nsoil-1 50 51 ! If we have permanent ice, we use the water ice thermal inertia from ground to last layer. 51 mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa 52 mthermdiff(ig,ik,islope)=inert_h2o_ice* 53 & inert_h2o_ice/volcapa 52 54 enddo 53 55 else 54 56 do ik=0,nsoil-1 55 mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa 57 mthermdiff(ig,ik,islope)=therm_i(ig,ik+1,islope)* 58 & therm_i(ig,ik+1,islope)/volcapa 56 59 enddo 57 60 endif 58 61 enddo 62 enddo 59 63 60 64 #ifdef MESOSCALE 61 65 do ig=1,ngrid 62 if ( therm_i(ig,1) .ge. inert_h2o_ice ) then 63 print *, "limit max TI ", therm_i(ig,1), inert_h2o_ice 66 do islope = 1,nslope 67 if ( therm_i(ig,1,islope) .ge. inert_h2o_ice ) then 68 print *, "limit max TI ", therm_i(ig,1,islope), inert_h2o_ice 64 69 do ik=0,nsoil-1 65 mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa 70 mthermdiff(ig,ik,islope)=inert_h2o_ice* 71 & inert_h2o_ice/volcapa 66 72 enddo 67 73 endif 74 enddo 68 75 enddo 69 76 #endif … … 71 78 ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities 72 79 do ig=1,ngrid 80 do islope = 1,nslope 73 81 do ik=1,nsoil-1 74 thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik) 75 & +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1)) 82 thermdiff(ig,ik,islope)=((layer(ik)-mlayer(ik-1)) 83 & *mthermdiff(ig,ik,islope) 84 & +(mlayer(ik)-layer(ik)) 85 & *mthermdiff(ig,ik-1,islope)) 76 86 & /(mlayer(ik)-mlayer(ik-1)) 77 87 ! write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik) 78 88 enddo 89 enddo 79 90 enddo 80 91 … … 92 103 93 104 do ig=1,ngrid 105 do islope = 1,nslope 94 106 ! d_k 95 107 do ik=1,nsoil-1 96 coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1)) 108 coefd(ig,ik,islope)=thermdiff(ig,ik,islope) 109 & /(mlayer(ik)-mlayer(ik-1)) 97 110 enddo 98 111 99 112 ! alph_{N-1} 100 alph(ig,nsoil-1 )=coefd(ig,nsoil-1)/101 & (coefq(nsoil-1)+coefd(ig,nsoil-1 ))113 alph(ig,nsoil-1,islope)=coefd(ig,nsoil-1,islope)/ 114 & (coefq(nsoil-1)+coefd(ig,nsoil-1,islope)) 102 115 ! alph_k 103 116 do ik=nsoil-2,1,-1 104 alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)* 105 & (1.-alph(ig,ik+1))+coefd(ig,ik)) 117 alph(ig,ik,islope)=coefd(ig,ik,islope)/ 118 & (coefq(ik)+coefd(ig,ik+1,islope)* 119 & (1.-alph(ig,ik+1,islope))+coefd(ig,ik,islope)) 106 120 enddo 107 121 108 122 ! capcal 109 123 ! Cstar 110 capcal(ig )=volcapa*layer(1)+111 & (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*112 & (timestep*(1.-alph(ig,1)))124 capcal(ig,islope)=volcapa*layer(1)+ 125 & (thermdiff(ig,1,islope)/(mlayer(1)-mlayer(0)))* 126 & (timestep*(1.-alph(ig,1,islope))) 113 127 ! Cs 114 capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* 115 & thermdiff(ig,1)/mthermdiff(ig,0)) 128 capcal(ig,islope)=capcal(ig,islope)/ 129 & (1.+mu*(1.0-alph(ig,1,islope))* 130 & thermdiff(ig,1,islope)/mthermdiff(ig,0,islope)) 116 131 ! write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig) 132 enddo ! islope 117 133 enddo ! of do ig=1,ngrid 118 134 … … 122 138 IF (.not.firstcall) THEN 123 139 ! First layer: 124 do ig=1,ngrid 125 tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)* 126 & thermdiff(ig,1)/mthermdiff(ig,0))/ 127 & (1.+mu*(1.0-alph(ig,1))* 128 & thermdiff(ig,1)/mthermdiff(ig,0)) 140 do islope = 1,nslope 141 do ig=1,ngrid 142 tsoil(ig,1,islope)=(tsurf(ig,islope)+mu*beta(ig,1,islope)* 143 & thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))/ 144 & (1.+mu*(1.0-alph(ig,1,islope))* 145 & thermdiff(ig,1,islope)/mthermdiff(ig,0,islope)) 129 146 enddo 130 147 ! Other layers: 131 148 do ik=1,nsoil-1 132 149 do ig=1,ngrid 133 tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik) 134 enddo 135 enddo 136 150 tsoil(ig,ik+1,islope)=alph(ig,ik,islope)* 151 & tsoil(ig,ik,islope)+beta(ig,ik,islope) 152 enddo 153 enddo 154 enddo ! islope 137 155 ENDIF! of if (.not.firstcall) 138 156 139 157 ! 2. Compute beta coefficients (preprocessing for next time step) 140 158 ! Bottom layer, beta_{N-1} 141 do ig=1,ngrid 142 beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) 143 & /(coefq(nsoil-1)+coefd(ig,nsoil-1)) 144 & +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1)) 159 do islope = 1,nslope 160 do ig=1,ngrid 161 beta(ig,nsoil-1,islope)=coefq(nsoil-1)*tsoil(ig,nsoil,islope) 162 & /(coefq(nsoil-1)+coefd(ig,nsoil-1,islope)) 163 & +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1,islope)) 145 164 enddo 146 165 … … 148 167 do ik=nsoil-2,1,-1 149 168 do ig=1,ngrid 150 beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+ 151 & coefd(ig,ik+1)*beta(ig,ik+1))/ 152 & (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) 153 & +coefd(ig,ik)) 169 beta(ig,ik,islope)=(coefq(ik)*tsoil(ig,ik+1,islope)+ 170 & coefd(ig,ik+1,islope)*beta(ig,ik+1,islope))/ 171 & (coefq(ik)+coefd(ig,ik+1,islope)* 172 & (1.0-alph(ig,ik+1,islope)) 173 & +coefd(ig,ik,islope)) 154 174 enddo 155 175 enddo … … 162 182 ! & (timestep*(1.-alph(ig,1))) 163 183 ! Fstar 164 fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* 165 & (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1)) 184 fluxgrd(ig,islope)=(thermdiff(ig,1,islope)/ 185 & (mlayer(1)-mlayer(0)))* 186 & (beta(ig,1,islope)+(alph(ig,1,islope)-1.0) 187 & *tsoil(ig,1,islope)) 166 188 167 189 ! mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0)) … … 169 191 ! & thermdiff(ig,1)/mthermdiff(ig,0)) 170 192 ! Fs 171 fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)* 172 & (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))* 173 & thermdiff(ig,1)/mthermdiff(ig,0)) 174 & -tsurf(ig)-mu*beta(ig,1)* 175 & thermdiff(ig,1)/mthermdiff(ig,0)) 176 enddo 177 193 fluxgrd(ig,islope)=fluxgrd(ig,islope)+(capcal(ig,islope) 194 & /timestep)* 195 & (tsoil(ig,1,islope)* 196 & (1.+mu*(1.0-alph(ig,1,islope))* 197 & thermdiff(ig,1,islope)/mthermdiff(ig,0,islope)) 198 & -tsurf(ig,islope)-mu*beta(ig,1,islope)* 199 & thermdiff(ig,1,islope)/mthermdiff(ig,0,islope)) 200 enddo 201 enddo ! islope 178 202 end 179 203 -
trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
r2826 r2900 44 44 45 45 !! variables 46 REAL,SAVE,ALLOCATABLE :: tsurf(: ) ! Surface temperature (K)47 REAL,SAVE,ALLOCATABLE :: emis(: ) ! Thermal IR surface emissivity48 REAL,SAVE,ALLOCATABLE :: capcal(: ) ! surface heat capacity (J m-2 K-1)49 REAL,SAVE,ALLOCATABLE :: fluxgrd(: ) ! surface conduction flux (W.m-2)50 REAL,ALLOCATABLE,SAVE :: qsurf(:,: ) ! tracer on surface (e.g. kg.m-2)51 REAL,SAVE,ALLOCATABLE :: watercap(: ) ! Surface water ice (kg.m-2)46 REAL,SAVE,ALLOCATABLE :: tsurf(:,:) ! Surface temperature (K) 47 REAL,SAVE,ALLOCATABLE :: emis(:,:) ! Thermal IR surface emissivity 48 REAL,SAVE,ALLOCATABLE :: capcal(:,:) ! surface heat capacity (J m-2 K-1) 49 REAL,SAVE,ALLOCATABLE :: fluxgrd(:,:) ! surface conduction flux (W.m-2) 50 REAL,ALLOCATABLE,SAVE :: qsurf(:,:,:) ! tracer on surface (e.g. kg.m-2) 51 REAL,SAVE,ALLOCATABLE :: watercap(:,:) ! Surface water ice (kg.m-2) 52 52 53 53 !$OMP THREADPRIVATE(tsurf,emis,capcal,fluxgrd,qsurf,watercap) … … 55 55 contains 56 56 57 subroutine ini_surfdat_h(ngrid,nq )57 subroutine ini_surfdat_h(ngrid,nq,nslope) 58 58 59 59 implicit none 60 60 integer,intent(in) :: ngrid ! number of atmospheric columns 61 61 integer,intent(in) :: nq ! number of tracers 62 62 integer,intent(in) :: nslope ! number of sub-grid scale slope 63 63 allocate(albedodat(ngrid)) 64 64 allocate(phisfi(ngrid)) … … 71 71 allocate(zthe(ngrid)) 72 72 allocate(z0(ngrid)) 73 allocate(qsurf(ngrid,nq ))74 allocate(tsurf(ngrid ))75 allocate(watercap(ngrid ))76 allocate(emis(ngrid ))77 allocate(capcal(ngrid ))78 allocate(fluxgrd(ngrid ))73 allocate(qsurf(ngrid,nq,nslope)) 74 allocate(tsurf(ngrid,nslope)) 75 allocate(watercap(ngrid,nslope)) 76 allocate(emis(ngrid,nslope)) 77 allocate(capcal(ngrid,nslope)) 78 allocate(fluxgrd(ngrid,nslope)) 79 79 allocate(hmons(ngrid)) 80 80 allocate(summit(ngrid)) -
trunk/LMDZ.MARS/libf/phymars/writediagfi.F
r2573 r2900 78 78 79 79 integer,save :: zitau=0 80 character(len=2 0),save :: firstnom='1234567890'80 character(len=27),save :: firstnom='1234567890' 81 81 !$OMP THREADPRIVATE(zitau,firstnom) 82 82
Note: See TracChangeset
for help on using the changeset viewer.
