Changeset 2935 for trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars
- Timestamp:
- Apr 12, 2023, 1:44:55 PM (21 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F
r2881 r2935 23 23 USE comvert_mod, ONLY: ap,bp,aps,bps,preff 24 24 USE comconst_mod, ONLY: kappa,g,pi 25 use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, 26 & subslope_dist,end_comslope_h,ini_comslope_h 27 use netcdf 25 28 implicit none 26 29 … … 44 47 REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1) 45 48 REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) 46 REAL,INTENT(OUT) :: tsurf(ngrid ) ! surface temperature47 REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx ) ! soil temperature48 REAL,INTENT(OUT) :: albedo(ngrid,2 ) ! surface albedo49 REAL,INTENT(OUT) :: emis(ngrid ) ! ground emissivity50 REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot )49 REAL,INTENT(OUT) :: tsurf(ngrid,nslope) ! surface temperature 50 REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx,nslope) ! soil temperature 51 REAL,INTENT(OUT) :: albedo(ngrid,2,nslope) ! surface albedo 52 REAL,INTENT(OUT) :: emis(ngrid,nslope) ! ground emissivity 53 REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot,nslope) 51 54 REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor 52 55 REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! sub grid cloud fraction 53 REAL,INTENT(OUT) :: watercap(ngrid ) ! infinite polar cap56 REAL,INTENT(OUT) :: watercap(ngrid,nslope) ! infinite polar cap 54 57 REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1) 55 58 REAL,INTENT(OUT) :: t(iip1,jjp1,llm) … … 128 131 c------------------------------------------------------ 129 132 real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) 130 real tsurfS(iip1,jjp1 ),tsoilS(iip1,jjp1,nsoilmx)133 real tsurfS(iip1,jjp1,nslope),tsoilS(iip1,jjp1,nsoilmx,nslope) 131 134 real inertiedatS(iip1,jjp1,nsoilmx) 132 real co2iceS(iip1,jjp1 ),emisS(iip1,jjp1)133 REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot )135 real co2iceS(iip1,jjp1,nslope),emisS(iip1,jjp1,nslope) 136 REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot,nslope) 134 137 real tauscalingS(iip1,jjp1) 135 138 real totcloudfracS(iip1,jjp1) 136 real watercapS(iip1,jjp1 )137 real albedoS(iip1,jjp1 )139 real watercapS(iip1,jjp1,nslope) 140 real albedoS(iip1,jjp1,nslope) 138 141 139 142 real ptotal, co2icetotal … … 153 156 real, dimension(:), allocatable :: mlayerold 154 157 real, dimension(:,:,:), allocatable :: uold,vold,told,q2old 155 real, dimension(:,:,:), allocatable :: tsoilold,qsurfold 156 real, dimension(:,:,:),allocatable :: tsoiloldnew 158 real, dimension(:,:,:,:), allocatable :: tsoilold,qsurfold 159 real, dimension(:,:,:),allocatable :: tsoiloldnew_noslope 160 real, dimension(:,:,:), allocatable :: tsoilold_noslope 161 real, dimension(:,:,:), allocatable ::qsurfold_noslope 162 real, dimension(:,:,:,:),allocatable :: tsoiloldnew 157 163 ! tsoiloldnew: old soil values, but along new subterranean grid 158 164 real, dimension(:,:,:), allocatable :: inertiedatold … … 160 166 real, dimension(:,:,:), allocatable :: inertiedatoldnew 161 167 real, dimension(:,:), allocatable :: psold,phisold 162 real, dimension(:,:), allocatable :: co2iceold,tsurfold 163 real, dimension(:,:), allocatable :: emisold 168 real, dimension(:,:,:), allocatable :: co2iceold,tsurfold 169 real, dimension(:,:), allocatable :: co2iceold_noslope 170 real, dimension(:,:), allocatable :: tsurfold_noslope 171 real, dimension(:,:,:), allocatable :: emisold 172 real, dimension(:,:), allocatable :: emisold_noslope 164 173 real, dimension(:,:,:,:), allocatable :: qold 165 174 real, dimension(:,:), allocatable :: tauscalingold 166 175 real, dimension(:,:), allocatable :: totcloudfracold 167 real, dimension(:,:), allocatable :: watercapold 168 real, dimension(:,:), allocatable :: albedoold 176 real, dimension(:,:,:), allocatable :: watercapold 177 real, dimension(:,:), allocatable :: watercapold_noslope 178 real, dimension(:,:,:), allocatable :: albedoold 179 real, dimension(:,:), allocatable :: albedoold_noslope 180 169 181 170 182 real tab_cntrl(100) … … 200 212 ! start_archive) and then .false. 201 213 logical :: old_co2ice=.false. 214 215 ! flag to check if nslope is not present (old start_archive) 216 integer :: nslope_read 217 logical :: no_slope=.false. 218 integer :: ndims 219 integer, dimension(:), allocatable :: dimids 202 220 c======================================================================= 203 221 … … 324 342 depthinterpol=.true. 325 343 endif 344 345 346 ! 1.2.1 find out the # of subgrid scale slope 347 ierr= NF_INQ_DIMID(nid,"nslope",dimid) 348 if (ierr.ne.NF_NOERR) then 349 write(*,*) "nslope not present, old startarchive" 350 write(*,*) "nslope=1" 351 nslope=1 352 no_slope=.true. 353 else 354 write(*,*) "nslope present" 355 ierr= NF_INQ_DIMLEN(nid,dimid,nslope_read) 356 if (ierr.ne.NF_NOERR) then 357 write(*,*) 'lect_start_archive error: cannot 358 & find nslope length' 359 stop 360 else 361 write(*,*) "lect_start_archive: nslope_read=",nslope_read 362 if(nslope_read.ne.1) then 363 write(*,*) 'lect_start_archive only works with nslope=1' 364 stop 365 else 366 nslope=1 367 no_slope=.false. 368 endif 369 endif 370 endif 371 def_slope(1)=-90. 372 def_slope(2)=90. 373 def_slope_mean(1)=0. 374 subslope_dist(:,:)=1. 326 375 327 376 ! 1.3 Report dimensions … … 331 380 write(*,*) "latitude: ",jmold 332 381 write(*,*) "altitude: ",lmold 382 write(*,*) "nslope: ",nslope 333 383 if (nqold.gt.0) then 334 384 write(*,*) "old tracers q*: ",nqold … … 358 408 allocate(phisold(imold+1,jmold+1)) 359 409 allocate(qold(imold+1,jmold+1,lmold,nqtot)) 360 allocate(co2iceold(imold+1,jmold+1 ))361 allocate(tsurfold(imold+1,jmold+1 ))362 allocate(emisold(imold+1,jmold+1 ))410 allocate(co2iceold(imold+1,jmold+1,nslope)) 411 allocate(tsurfold(imold+1,jmold+1,nslope)) 412 allocate(emisold(imold+1,jmold+1,nslope)) 363 413 allocate(q2old(imold+1,jmold+1,lmold+1)) 364 414 ! allocate(tsoilold(imold+1,jmold+1,nsoilmx)) 365 allocate(tsoilold(imold+1,jmold+1,nsoilold ))366 allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx ))415 allocate(tsoilold(imold+1,jmold+1,nsoilold,nslope)) 416 allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx,nslope)) 367 417 allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia 368 418 allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx)) … … 370 420 allocate(surfithold(imold+1,jmold+1)) 371 421 allocate(mlayerold(nsoilold)) 372 allocate(qsurfold(imold+1,jmold+1,nqtot ))422 allocate(qsurfold(imold+1,jmold+1,nqtot,nslope)) 373 423 allocate(tauscalingold(imold+1,jmold+1)) 374 424 allocate(totcloudfracold(imold+1,jmold+1)) 375 allocate(watercapold(imold+1,jmold+1)) 376 allocate(albedoold(imold+1,jmold+1)) 425 allocate(watercapold(imold+1,jmold+1,nslope)) 426 allocate(albedoold(imold+1,jmold+1,nslope)) 427 428 if(no_slope) then 429 allocate(co2iceold_noslope(imold+1,jmold+1)) 430 allocate(tsurfold_noslope(imold+1,jmold+1)) 431 allocate(emisold_noslope(imold+1,jmold+1)) 432 allocate(watercapold_noslope(imold+1,jmold+1)) 433 allocate(albedoold_noslope(imold+1,jmold+1)) 434 allocate(qsurfold_noslope(imold+1,jmold+1,nqtot)) 435 allocate(tsoilold_noslope(imold+1,jmold+1,nsoilold)) 436 allocate(tsoiloldnew_noslope(imold+1,jmold+1,nsoilmx)) 437 endif 377 438 378 439 allocate(var (imold+1,jmold+1,llm)) … … 722 783 ENDIF 723 784 ENDIF 724 #ifdef NC_DOUBLE 725 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold) 726 #else 727 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold) 728 #endif 785 if(no_slope) then 786 ierr = nf90_get_var(nid, nvarid,co2iceold_noslope) 787 co2iceold(:,:,1)=co2iceold_noslope(:,:) 788 else 789 ierr = nf90_get_var(nid, nvarid,co2iceold) 790 endif 729 791 IF (ierr .NE. NF_NOERR) THEN 730 792 PRINT*, "lect_start_archive: Failed loading <co2ice>" … … 738 800 CALL abort 739 801 ENDIF 740 #ifdef NC_DOUBLE 741 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold) 742 #else 743 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold) 744 #endif 802 if(no_slope) then 803 ierr = nf90_get_var(nid, nvarid,emisold_noslope) 804 emisold(:,:,1)=emisold_noslope(:,:) 805 else 806 ierr = nf90_get_var(nid, nvarid,emisold) 807 endif 745 808 IF (ierr .NE. NF_NOERR) THEN 746 809 PRINT*, "lect_start_archive: Failed loading <emis>" … … 768 831 CALL abort 769 832 ENDIF 770 #ifdef NC_DOUBLE 771 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold) 772 #else 773 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold) 774 #endif 833 if(no_slope) then 834 ierr = nf90_get_var(nid, nvarid,tsurfold_noslope) 835 tsurfold(:,:,1)=tsurfold_noslope(:,:) 836 else 837 ierr = nf90_get_var(nid, nvarid,tsurfold) 838 endif 775 839 IF (ierr .NE. NF_NOERR) THEN 776 840 PRINT*, "lect_start_archive: Failed loading <tsurf>" … … 832 896 IF (ierr .NE. NF_NOERR) THEN 833 897 PRINT*, "lect_start_archive: <watercap> not in file" 834 watercapold(:,: ) = 0.898 watercapold(:,:,:) = 0. 835 899 ELSE 836 #ifdef NC_DOUBLE 837 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count, 838 & watercapold) 839 #else 840 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count, 841 & watercapold) 842 #endif 900 if(no_slope) then 901 ierr = nf90_get_var(nid, nvarid,watercapold_noslope) 902 watercapold(:,:,1)=watercapold_noslope(:,:) 903 else 904 ierr = nf90_get_var(nid, nvarid,watercapold) 905 endif 843 906 IF (ierr .NE. NF_NOERR) THEN 844 907 PRINT*, "lect_start_archive: Failed loading <watercap>" … … 851 914 IF (ierr .NE. NF_NOERR) THEN 852 915 PRINT*, "lect_start_archive: <albedo> not in file" 853 albedoold(:,: ) = 0.916 albedoold(:,:,:) = 0. 854 917 ELSE 855 #ifdef NC_DOUBLE 856 ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,count, 857 & albedoold) 858 #else 859 ierr = NF_GET_VARA_REAL(nid,nvarid,start,count, 860 & albedoold) 861 #endif 918 if(no_slope) then 919 ierr = nf90_get_var(nid, nvarid,albedoold_noslope) 920 albedoold(:,:,1)=albedoold_noslope(:,:) 921 else 922 ierr = nf90_get_var(nid, nvarid,albedoold) 923 endif 862 924 IF (ierr .NE. NF_NOERR) THEN 863 925 PRINT*, "lect_start_archive: Failed loading <albedo>" … … 876 938 c ------------------------------------------- 877 939 ! Surface tracers: 878 qsurfold(1:imold+1,1:jmold+1,1:nqtot )=0940 qsurfold(1:imold+1,1:jmold+1,1:nqtot,1:nslope)=0 879 941 880 942 DO iq=1,nqtot … … 906 968 print*, "which (constant) value should it be initialized to?" 907 969 read(*,*) tmpval 908 qsurfold(1:imold+1,1:jmold+1,iq )=tmpval970 qsurfold(1:imold+1,1:jmold+1,iq,1:nslope)=tmpval 909 971 ELSE ! tracer exists in file, load it 910 #ifdef NC_DOUBLE 911 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count, 912 & qsurfold(1,1,iq)) 913 #else 914 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count, 915 & qsurfold(1,1,iq)) 916 #endif 972 if(no_slope) then 973 ierr = nf90_get_var(nid, nvarid,qsurfold_noslope(:,:,iq)) 974 qsurfold(:,:,iq,1)=qsurfold_noslope(:,:,iq) 975 else 976 ierr = nf90_get_var(nid, nvarid,qsurfold(:,:,iq,:)) 977 endif 917 978 IF (ierr .NE. NF_NOERR) THEN 918 979 PRINT*, "lect_start_archive: ", … … 946 1007 CALL abort 947 1008 ENDIF 948 #ifdef NC_DOUBLE 949 ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count, 950 & tsoilold(1,1,isoil)) 951 #else 952 ierr = NF_GET_VARA_REAL(nid, nvarid,start,count, 953 & tsoilold(1,1,isoil)) 954 #endif 1009 if(no_slope) then 1010 ierr = nf90_get_var(nid, nvarid,tsoilold_noslope(:,:,isoil)) 1011 tsoilold(:,:,:,1)=tsoilold_noslope(:,:,:) 1012 else 1013 ierr = nf90_get_var(nid, nvarid,tsoilold(:,:,isoil,:)) 1014 endif 955 1015 IF (ierr .NE. NF_NOERR) THEN 956 1016 PRINT*, "lect_start_archive: ", … … 972 1032 call abort 973 1033 else 974 #ifdef NC_DOUBLE 975 ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold) 976 #else 977 ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold) 978 #endif 1034 if(no_slope) then 1035 ierr = nf90_get_var(nid, nvarid,tsoilold_noslope) 1036 tsoilold(:,:,:,1)=tsoilold_noslope(:,:,:) 1037 else 1038 ierr = nf90_get_var(nid, nvarid,tsoilold) 1039 endif 979 1040 endif ! of if (ierr.ne.NF_NOERR) 980 1041 … … 1130 1191 1131 1192 c Glace CO2 1132 call interp_horiz (co2iceold ,co2ices,imold,jmold,iim,jjm,1,1133 & 1193 call interp_horiz (co2iceold(:,:,1),co2ices(:,:,1), 1194 & imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1134 1195 1135 1196 c Temperature de surface 1136 call interp_horiz (tsurfold ,tsurfs,imold,jmold,iim,jjm,1,1137 & 1138 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs ,tsurf)1197 call interp_horiz (tsurfold(:,:,1),tsurfs(:,:,1) 1198 & ,imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1199 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs(:,:,1),tsurf(:,1)) 1139 1200 c write(44,*) 'tsurf', tsurf 1140 1201 … … 1147 1208 1148 1209 c Emissivite de la surface 1149 call interp_horiz (emisold ,emiss,imold,jmold,iim,jjm,1,1150 & 1151 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss ,emis)1210 call interp_horiz (emisold(:,:,1),emiss(:,:,1) 1211 & ,imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1212 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss(:,:,1),emis(:,1)) 1152 1213 1153 1214 c Dust conversion factor … … 1162 1223 1163 1224 c Watercap 1164 call interp_horiz (watercapold,watercaps,imold,jmold,iim, 1165 & jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1166 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,watercaps,watercap) 1225 call interp_horiz (watercapold(:,:,1),watercaps(:,:,1), 1226 & imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1227 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,watercaps(:,:,1), 1228 & watercap(:,1)) 1167 1229 1168 1230 c Surface albedo 1169 call interp_horiz (albedoold ,albedoS,imold,jmold,iim,1170 & 1171 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,albedoS ,albedo(:,1))1231 call interp_horiz (albedoold(:,:,1),albedoS(:,:,1), 1232 & imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv) 1233 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,albedoS(:,:,1),albedo(:,1,1)) 1172 1234 1173 albedo(:,2 )=albedo(:,1)1235 albedo(:,2,1)=albedo(:,1,1) 1174 1236 1175 1237 c write(46,*) 'emis',emis … … 1191 1253 DO i=1,iim 1192 1254 ptotal=ptotal+ps(i,j)*aire(i,j)/g 1193 co2icetotal = co2icetotal + co2iceS(i,j )*aire(i,j)1255 co2icetotal = co2icetotal + co2iceS(i,j,1)*aire(i,j) 1194 1256 END DO 1195 1257 END DO … … 1217 1279 DO j=1,jjp1 1218 1280 DO i=1,iip1 1219 co2iceS(i,j )=co2iceS(i,j) *co2icetotalold/co2icetotal1281 co2iceS(i,j,1)=co2iceS(i,j,1)*co2icetotalold/co2icetotal 1220 1282 END DO 1221 1283 END DO … … 1304 1366 do j=1,jmold+1 1305 1367 ! copy values 1306 oldval(1)=tsurfold(i,j )1307 oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold )1368 oldval(1)=tsurfold(i,j,1) 1369 oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold,1) 1308 1370 ! build vertical coordinate 1309 1371 oldgrid(1)=0. ! ground … … 1317 1379 & mlayer,newval,nsoilmx) 1318 1380 ! copy result in tsoilold 1319 tsoiloldnew(i,j,1:nsoilmx )=newval(1:nsoilmx)1381 tsoiloldnew(i,j,1:nsoilmx,1)=newval(1:nsoilmx) 1320 1382 enddo 1321 1383 enddo … … 1336 1398 do j=1,jmold+1 1337 1399 ! copy values 1338 oldval(1)=tsurfold(i,j )1339 oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold )1400 oldval(1)=tsurfold(i,j,1) 1401 oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold,1) 1340 1402 ! interpolate 1341 1403 call interp_line(oldgrid,oldval,nsoilold+1, 1342 1404 & mlayer,newval,nsoilmx) 1343 1405 ! copy result in tsoilold 1344 tsoiloldnew(i,j,1:nsoilmx )=newval(1:nsoilmx)1406 tsoiloldnew(i,j,1:nsoilmx,1)=newval(1:nsoilmx) 1345 1407 enddo 1346 1408 enddo … … 1352 1414 1353 1415 else 1354 tsoiloldnew(:,:,: )=tsoilold(:,:,:)1416 tsoiloldnew(:,:,:,1)=tsoilold(:,:,:,1) 1355 1417 endif ! of if (depthinterpol) 1356 1418 endif ! of if (olddepthdef) 1357 1419 1358 1420 ! Do the horizontal interpolation 1359 call interp_horiz(tsoiloldnew ,tsoilS,1421 call interp_horiz(tsoiloldnew(:,:,:,1),tsoilS(:,:,:,1), 1360 1422 & imold,jmold,iim,jjm,nsoilmx, 1361 1423 & rlonuold,rlatvold,rlonu,rlatv) 1362 1424 1363 1425 ! Reshape tsoilS to scalar grid as tsoil 1364 call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil) 1426 call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS(:,:,:,1), 1427 & tsoil(:,:,:)) 1365 1428 1366 1429 … … 1432 1495 c traceurs surface 1433 1496 do iq = 1, nqtot 1434 call interp_horiz(qsurfold(1,1,iq ) ,qsurfs(1,1,iq),1497 call interp_horiz(qsurfold(1,1,iq,1) ,qsurfs(1,1,iq,1), 1435 1498 & imold,jmold,iim,jjm,1, 1436 1499 & rlonuold,rlatvold,rlonu,rlatv) 1437 1500 enddo 1438 1501 1439 call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf) 1502 call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs(:,:,:,1), 1503 & qsurf(:,:,1)) 1440 1504 1441 1505 c traceurs 3D … … 1488 1552 enddo 1489 1553 1490 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,qsurf(:,igcm_co2)) 1554 call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices(:,:,1), 1555 & qsurf(:,igcm_co2,1)) 1491 1556 1492 1557 c-----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.