- Timestamp:
- Oct 15, 2020, 11:46:56 AM (4 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 2 edited
-
README (modified) (1 diff)
-
util/solzenangle.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2415 r2416 3164 3164 "tau_pref_gcm" and "tau_pref_scenario" which are repectively the visible 3165 3165 dust opacity column (at 610 Pa) in the GCM and prescribed by a dust scenario. 3166 3167 == 15/10/2020 == AB 3168 Correction in solzenangle.F90 : the missing value was not initialized when the attribute was not 3169 present in the input file. 3170 + Resolved Ticket #61 : add comments in the code and in the output file to make the program more 3171 transparent for the users -
trunk/LMDZ.MARS/util/solzenangle.F90
r2289 r2416 26 26 ! filename(): output file name 27 27 ! vartmp(): temporary variable name (read from netcdf input file) 28 !character (len=1) :: ccopy 29 ! ccpy: 'y' or 'n' answer28 character (len=500) :: globatt 29 ! globatt: [netcdf] global attribute of the output file 30 30 31 31 character (len=50) :: altlong_name,altunits,altpositive … … 100 100 ! var3D_lt(,,,): 4D array to store a field in local time coordinate 101 101 real, dimension(:), allocatable :: lt_gcm 102 real, dimension(:), allocatable :: lt_outc 102 real, dimension(:), allocatable :: lt_outc ! local lt_out, put in sol decimal value 103 103 real, dimension(:), allocatable :: var_gcm 104 104 real, dimension(:), allocatable :: intsol 105 105 real, dimension(:,:,:), allocatable :: lt_out 106 ! lt_out: computed local time (hour) corresponding to sza, put in output file 106 107 real :: sza ! solar zenith angle 107 108 character (len=30) :: szastr ! to store the sza value in a string … … 237 238 !============================================================================== 238 239 239 240 241 242 243 244 245 246 240 write(*,*) 241 write(*,*) "opening "//trim(file)//"..." 242 ierr = NF_OPEN(file,NF_NOWRITE,nid) 243 if (ierr.NE.NF_NOERR) then 244 write(*,*) 'ERROR: Pb opening file '//trim(file) 245 write(*,*) NF_STRERROR(ierr) 246 stop "" 247 endif 247 248 248 249 !============================================================================== … … 311 312 312 313 ! First call; initialize/allocate 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 #ifdef NC_DOUBLE 335 336 #else 337 338 #endif 339 340 341 342 343 344 345 #ifdef NC_DOUBLE 346 347 #else 348 349 #endif 350 351 352 353 354 355 356 #ifdef NC_DOUBLE 357 358 #else 359 360 #endif 361 362 363 364 365 314 allocate(lat(latlen),stat=ierr) 315 if (ierr.ne.0) then 316 write(*,*) "Error: failed to allocate lat(latlen)" 317 stop 318 endif 319 allocate(lon(lonlen),stat=ierr) 320 if (ierr.ne.0) then 321 write(*,*) "Error: failed to allocate lon(lonlen)" 322 stop 323 endif 324 allocate(alt(altlen),stat=ierr) 325 if (ierr.ne.0) then 326 write(*,*) "Error: failed to allocate alt(altlen)" 327 stop 328 endif 329 allocate(ctl(ctllen),stat=ierr) 330 if (ierr.ne.0) then 331 write(*,*) "Error: failed to allocate ctl(ctllen)" 332 stop 333 endif 334 335 #ifdef NC_DOUBLE 336 ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat) 337 #else 338 ierr = NF_GET_VAR_REAL(nid,latvar,lat) 339 #endif 340 if (ierr.ne.0) then 341 write(*,*) "Error: failed to load latitude" 342 write(*,*) NF_STRERROR(ierr) 343 stop 344 endif 345 346 #ifdef NC_DOUBLE 347 ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon) 348 #else 349 ierr = NF_GET_VAR_REAL(nid,lonvar,lon) 350 #endif 351 if (ierr.ne.0) then 352 write(*,*) "Error: failed to load longitude" 353 write(*,*) NF_STRERROR(ierr) 354 stop 355 endif 356 357 #ifdef NC_DOUBLE 358 ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt) 359 #else 360 ierr = NF_GET_VAR_REAL(nid,altvar,alt) 361 #endif 362 if (ierr.ne.0) then 363 write(*,*) "Error: failed to load altitude" 364 write(*,*) NF_STRERROR(ierr) 365 stop 366 endif 366 367 ! Get altitude attributes to handle files with any altitude type 367 ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)368 ierr=nf_get_att_text(nid,altvar,'units',altunits)369 ierr=nf_get_att_text(nid,altvar,'positive',altpositive)370 371 372 #ifdef NC_DOUBLE 373 374 #else 375 376 #endif 377 378 379 380 381 382 368 ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name) 369 ierr=nf_get_att_text(nid,altvar,'units',altunits) 370 ierr=nf_get_att_text(nid,altvar,'positive',altpositive) 371 372 if (ctllen .ne. 0) then 373 #ifdef NC_DOUBLE 374 ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl) 375 #else 376 ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl) 377 #endif 378 if (ierr.ne.0) then 379 write(*,*) "Error: failed to load controle" 380 write(*,*) NF_STRERROR(ierr) 381 stop 382 endif 383 endif ! of if (ctllen .ne. 0) 383 384 384 385 !============================================================================== … … 430 431 endif 431 432 432 433 433 !============================================================================== 434 434 ! 2.4.1 Select local times to be saved "Time" in output file 435 435 !============================================================================== 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 436 write(*,*) "Planet side of the Solar Zenith Angle ?" 437 read(*,*) planetside 438 if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then 439 write(*,*) "planetside = "//planetside 440 write(*,*) "Error: the planet side must be morning or evening" 441 stop 442 endif 443 write(*,*) "Solar Zenith angle to interpolate ?" 444 write(*,*) "(between 0 and 180 deg, 90deg for the terminator at the surface)" 445 read(*,*) sza 446 if ((sza.lt.0).or.(sza.ge.180)) then 447 write(*,*) "ERROR: the solar zenith angle must lie within [0;180[ deg" 448 stop 449 endif 450 write(*,*) "Beginning date of the simulation file " 451 write(*,*) "(or median date of the simulation for a stats)?" 452 write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)" 453 read(*,*) starttimeoffset 454 455 if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then 456 write(*,*) 'We detect a ""stats"" file' 457 stats=.true. 458 458 ! We rewrite the time from "stats" from 0 to 1 sol... 459 do it=1,timelen ! loop temps in460 461 462 463 464 459 do it=1,timelen ! loop time input file 460 time(it) = time(it)/24. 461 end do 462 nsol =1 463 else ! case of a diagfi or concatnc file 464 stats=.false. 465 465 ! Number of sol in input file 466 nsol = int(time(timelen)) - int(time(1)) 467 end if 468 469 allocate(intsol(nsol),stat=ierr) 470 if (ierr.ne.0) then 471 write(*,*) "Error: failed to allocate intsol(nsol)" 472 stop 473 endif 474 allocate(lt_out(lonlen,latlen,nsol),stat=ierr) 475 if (ierr.ne.0) then 476 write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)" 477 stop 478 endif 479 do it=1,nsol 480 intsol(it)= int(time(1)) + it-1. 481 do ilon=1,lonlen 482 ! For the computation of Ls, it's supposed that the morning/evening 466 nsol = int(time(timelen)) - int(time(1)) 467 end if 468 469 allocate(intsol(nsol),stat=ierr) 470 if (ierr.ne.0) then 471 write(*,*) "Error: failed to allocate intsol(nsol)" 472 stop 473 endif 474 allocate(lt_out(lonlen,latlen,nsol),stat=ierr) 475 if (ierr.ne.0) then 476 write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)" 477 stop 478 endif 479 do it=1,nsol 480 intsol(it)= int(time(1)) + it-1. 481 do ilon=1,lonlen 482 ! For the computation of Ls, we try to take into account the rotation time 483 ! of Mars it's supposed that the morning/evening 483 484 ! correspond to LT=6h/18h at the given longitude, which is 484 485 ! then reported to a sol value at longitude 0 485 486 487 488 489 490 491 486 if (trim(planetside).eq."morning") then 487 tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset 488 else !if trim(planetside).eq."evening" 489 tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset 490 endif 491 call sol2ls(tmpsol,tmpLs) 492 do ilat=1,latlen 492 493 ! Compute the Local Time of the solar zenith angle at a given longitude 493 494 call LTsolzenangle(lat(ilat),tmpLs,planetside,sza,miss_lt,tmpLT) 494 495 ! Deduce the sol value at this longitude 495 496 497 498 499 500 501 502 503 504 505 506 507 496 if (tmpLT.eq.miss_lt) then 497 ! if there is no LT corresponding to the sza, put a missing value in lt_out 498 lt_out(ilon,ilat,it)=miss_lt 499 else 500 lt_out(ilon,ilat,it)=tmpLT 501 endif 502 enddo 503 enddo 504 enddo 505 506 if (nsol.gt.1) then 507 timelen_tot=timelen 508 else 508 509 ! if only 1 sol available, we must add 1 timestep for the interpolation 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 510 timelen_tot=timelen+1 511 endif 512 allocate(lt_gcm(timelen_tot),stat=ierr) 513 if (ierr.ne.0) then 514 write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)" 515 stop 516 endif 517 allocate(var_gcm(timelen_tot),stat=ierr) 518 if (ierr.ne.0) then 519 write(*,*) "Error: failed to allocate var_gcm(timelen_tot)" 520 stop 521 endif 522 allocate(lt_outc(nsol),stat=ierr) 523 if (ierr.ne.0) then 524 write(*,*) "Error: failed to allocate lt_outc(nsol)" 525 stop 526 endif 526 527 527 528 !============================================================================== 528 529 ! 2.4.2. Get output file name 529 530 !============================================================================== 530 531 532 533 534 531 if (trim(planetside).eq."morning") then 532 filename=file(1:len_trim(file)-3)//"_MO.nc" 533 else ! if trim(planetside).eq."evening" 534 filename=file(1:len_trim(file)-3)//"_EV.nc" 535 endif 535 536 536 537 !============================================================================== … … 539 540 540 541 541 542 543 544 545 546 547 548 542 ! Initialize output file's lat,lon,alt and time dimensions 543 call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,& 544 altlong_name,altunits,altpositive,& 545 nout,latdimout,londimout,altdimout,timedimout,timevarout) 546 547 ! Initialize output file's aps,bps and phisinit variables 548 call init2(nid,lonlen,latlen,altlen,altdim,& 549 nout,londimout,latdimout,altdimout) 549 550 550 551 … … 640 641 endif 641 642 642 643 644 645 646 643 units=" " 644 title=" " 645 ierr=nf_get_att_text(nid,varid,"title",title) 646 ierr=nf_get_att_text(nid,varid,"units",units) 647 call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr) 647 648 648 649 … … 662 663 #endif 663 664 665 ! If there is no missing value attribute in the input variable, create one 666 if (miss.ne.NF_NOERR) then 667 missing = miss_lt 668 miss = NF_NOERR 669 endif 670 664 671 !============================================================================== 665 672 ! 2.5.4 Interpolation in local time 666 673 !============================================================================== 667 674 668 do ilon=1,lonlen 669 ! write(*,*) 'processing longitude =', lon(ilon) 670 ! Local time at each longitude : 671 do it=1,timelen ! loop temps in 672 lt_gcm(it) = time(it) + lon(ilon)/360. 673 enddo 674 if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1 675 676 do ilat=1,latlen 675 do ilon=1,lonlen 676 ! Recast GCM local time (sol decimal value) at each longitude : 677 do it=1,timelen ! loop time input file 678 lt_gcm(it) = time(it) + lon(ilon)/360. 679 enddo 680 if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1 681 682 do ilat=1,latlen 683 do it=1,nsol ! loop time local time 684 if (lt_out(ilon,ilat,it).eq.miss_lt) then 685 lt_outc(it)=miss_lt 686 else 687 ! Put computed sza local time (hour) into a sol decimal value : 688 lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24. 689 endif 690 ! If the computed local time exceeds the GCM input file time bounds, 691 ! put a missing value in the local time variable 692 if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then 693 lt_out(ilon,ilat,it)=miss_lt 694 lt_outc(it)=miss_lt 695 endif 696 enddo ! local time 697 698 if (ndim==3) then 699 do it=1,timelen ! loop time input file 700 var_gcm(it) = var3d(ilon,ilat,it,1) 701 enddo 702 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1) 677 703 do it=1,nsol ! loop time local time 678 if (lt_out (ilon,ilat,it).eq.miss_lt) then679 lt_outc(it)=miss_lt704 if (lt_outc(it).eq.miss_lt) then 705 var3d_lt(ilon,ilat,it,1)=missing 680 706 else 681 lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24. 682 endif 683 ! If the computed local time exceeds the GCM input file time bounds, 684 ! put a missing value in the local time variable 685 if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then 686 lt_out(ilon,ilat,it)=miss_lt 687 lt_outc(it)=miss_lt 707 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),& 708 missing,lt_gcm,var_gcm,timelen_tot) 688 709 endif 689 710 enddo ! local time 690 711 691 if (ndim==3) then 692 do it=1,timelen ! loop temps in 693 var_gcm(it) = var3d(ilon,ilat,it,1) 712 else if (ndim==4) then 713 do ialt=1,altlen 714 do it=1,timelen ! loop time input file 715 var_gcm(it) = var3d(ilon,ilat,ialt,it) 694 716 enddo 695 717 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1) 696 718 do it=1,nsol ! loop time local time 697 719 if (lt_outc(it).eq.miss_lt) then 698 var3d_lt(ilon,ilat,i t,1)=missing720 var3d_lt(ilon,ilat,ialt,it)=missing 699 721 else 700 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,i t,1),&722 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),& 701 723 missing,lt_gcm,var_gcm,timelen_tot) 702 724 endif 703 enddo ! local time 704 705 else if (ndim==4) then 706 do ialt=1,altlen 707 do it=1,timelen ! loop temps in 708 var_gcm(it) = var3d(ilon,ilat,ialt,it) 709 enddo 710 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1) 711 do it=1,nsol ! loop time local time 712 if (lt_outc(it).eq.miss_lt) then 713 var3d_lt(ilon,ilat,ialt,it)=missing 714 else 715 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),& 716 missing,lt_gcm,var_gcm,timelen_tot) 717 endif 718 enddo ! local time 719 enddo ! alt 720 endif 721 722 enddo ! lat 723 end do ! lon 725 enddo ! local time 726 enddo ! alt 727 endif 728 729 enddo ! lat 730 end do ! lon 724 731 725 732 … … 739 746 endif 740 747 741 ! In case there is a "valid_range" attribute 742 ! Write "valid_range" and "missing_value" attributes in output file 743 if (miss.eq.NF_NOERR) then 744 call missing_value(nout,varidout,validr,miss,valid_range,missing) 745 endif 748 ! Write "missing_value" (and "valid_range" if there is one) attributes in output file 749 call missing_value(nout,varidout,validr,miss,valid_range,missing) 746 750 747 751 ! free var3d() array … … 809 813 ierr=nf_close(nid) 810 814 815 ! Write output file's global attribute 816 globatt = "GCM simulation that has been interpolated at the prescribed "//trim(planetside)//& 817 "-side solar zenith angle of "//trim(adjustl(szastr))//"deg all around the globe. "//& 818 "Note that the variable LTsolzenangle stands for the Local (at lon,lat) True Solar Time corresponding "//& 819 "to the SZA. You can derive the corresponding sol at lon=0deg by combining LTsolzenangle, longitude and "//& 820 "Time (integer sol value at lon=0deg) variables. For any question on the program, please contact "//& 821 "antoine.bierjon@lmd.jussieu.fr" 822 ! Switch to netcdf define mode 823 ierr = NF_REDEF (nout) 824 ierr = NF_PUT_ATT_TEXT(nout,NF_GLOBAL,'comment',len_trim(globatt),trim(globatt)) 825 ! End netcdf define mode 826 ierr = NF_ENDDEF(nout) 811 827 ! Close output file 812 828 ierr=NF_CLOSE(nout) … … 941 957 call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,& 942 958 (/timedimout/),timevarout,ierr) 959 ! Switch to netcdf define mode 960 ierr = NF_REDEF (nout) 961 ierr = NF_PUT_ATT_TEXT(nout,timevarout,'comment',135,& 962 "The true unit of the variable Time is the integer value of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret") 963 ! End netcdf define mode 964 ierr = NF_ENDDEF(nout) 943 965 944 966 !==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
