Ignore:
Timestamp:
Aug 22, 2016, 11:03:18 AM (8 years ago)
Author:
aslmd
Message:

Goal: make the mesoscale model compliant
with new interface philosophy and adapted
to easy interfacing with any physics

Commit 1: moved all operations specific to
one planet (Mars so far) in dedicated routines
-- update_inputs_physiq_mod.F: from WRF to LMD physics
-- update_outputs_physiq_mod.F: from LMD physics WRF
Those new routines are in the WRFV2/phys/ folder
Changed Makefile in this folder to compile those

Interfacing for a new planet (e.g. Venus)
will just require to adapt those files
to the new planet. This will be done later.

NB: mesoscale model still not working
with HEAD version
but changes in this commit verified with v1520 (compile+debug+run)

Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/Makefile

    r94 r1577  
    99        module_physics_addtendc.o \
    1010        module_physics_init.o \
     11        update_inputs_physiq_mod.o \
     12        update_outputs_physiq_mod.o \
    1113        module_lmd_driver.o
    1214 

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F.new

    r1411 r1577  
    5858   USE module_model_constants
    5959   USE module_wrf_error
    60 
    6160#ifndef NOPHYS
    62    !! MODULES FROM LMD PHYSICS
    63    use slope_mod, only: theta_sl, psi_sl
    64    use surfdat_h, only: phisfi, albedodat,  &
    65                         zmea, zstd, zsig, zgam, zthe, z0, &
    66                         emissiv,albedice,iceradius, &
    67                         emisice,dtemisice, &
    68                         z0_default, &
    69                         tsurf, co2ice, emis, qsurf
    70    use comsoil_h, only: inertiedat,mlayer,layer, &
    71                         volcapa, tsoil
    72    use planete_h, only: year_day,periheli,aphelie, &
    73                         peri_day,obliquit,emin_turb, &
    74                         lmixmin
    75    use turb_mod, only: q2,wstar,ustar,sensibFlux,&
    76                        hfmax_th,zmax_th,turb_resolved
    77    use dimradmars_mod, only: fluxrad
    78    use tracer_mod, only: noms
    79    use comcstfi_h, only: omeg,mugaz
    80    use comgeomfi_h, only: ini_fillgeom
    81    USE comm_wrf
     61   USE update_inputs_physiq_mod !! to set inputs for physiq
     62   USE update_outputs_physiq_mod !! to get outputs from physiq
     63   USE comm_wrf !! to get fields to be written from physiq
    8264#endif
    8365
    … …  
    146128   INTEGER ::    i,j,k,its,ite,jts,jte,ij,kk
    147129   INTEGER ::    subs,iii
    148    REAL :: lay1,alpha
    149130   
    150131
    … …  
    157138   INTEGER :: ngrid,nlayer,nq,nsoil
    158139   REAL :: pday,ptime,MY 
    159    REAL :: aire_val,lat_val,lon_val
    160    REAL :: phisfi_val,albedodat_val,inertiedat_val
    161    REAL :: zmea_val,zstd_val,zsig_val,zgam_val,zthe_val
    162    REAL :: theta_val, psi_val
     140   REAL :: phisfi_val
    163141   LOGICAL :: firstcall,lastcall,tracerdyn
    164    REAL,DIMENSION(:),ALLOCATABLE :: plon, plat, parea
    165    REAL,DIMENSION(:),ALLOCATABLE :: q2_val, qsurf_val, tsoil_val
    166    REAL :: wstar_val,fluxrad_val
    167    REAL,DIMENSION(:),ALLOCATABLE :: isoil_val, dsoil_val
    168    REAL :: z0_val
    169142   ! ----------
    170143   REAL,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,pw
    … …  
    213186!!!IDEALIZED IDEALIZED
    214187      REAL :: lat_input, lon_input, ls_input, lct_input
     188      LOGICAL :: isles
    215189!!!IDEALIZED IDEALIZED
    216190
    … …  
    225199!==================================================================
    226200
     201!! determine here if this is turbulence-resolving or not
    227202IF (JULYR .ne. 9999) THEN
    228     turb_resolved = .false.  ! "True" LES is not available in this version
     203    isles = .false.  ! "True" LES is not available in this version
    229204    PRINT *, '*** REAL-CASE SIMULATION ***'
    230205ELSE
    … …  
    234209          PRINT *, '*** diff_opt = 2 *** km_opt = 2'
    235210          PRINT *, '*** forcing is isfflx = ',isfflx
    236           turb_resolved = .true.
     211          isles = .true.
    237212          !! SPECIAL LES
    238213     ELSE
    … …  
    240215          PRINT *, '*** diff_opt = ',diff_opt
    241216          PRINT *, '*** km_opt = ',km_opt
    242           turb_resolved = .false.
     217          isles = .false.
    243218          !! IDEALIZED, no LES
    244219          !! cependant, ne veut-on pas pouvoir
    … …  
    258233jte = j_end(num_tiles)
    259234!!
    260 IF (turb_resolved .eqv. .false.) THEN
     235IF (isles .eqv. .false.) THEN
    261236 relax=0
    262237 sponge_top=0               ! another value than 0 triggers instabilities 
    … …  
    267242jps=jts
    268243jpe=jte
    269 IF (turb_resolved .eqv. .false.) THEN
     244IF (isles .eqv. .false.) THEN
    270245 IF (ips .eq. ids)   ips=its+relax !! IF tests necesary for parallel runs
    271246 IF (ipe .eq. ide-1) ipe=ite-relax
    … …  
    274249ENDIF
    275250kps=kts         !! start at surface
    276 IF (turb_resolved .eqv. .false.) THEN
     251IF (isles .eqv. .false.) THEN
    277252 kpe=kte-sponge_top
    278253ELSE
    … …  
    375350ENDIF
    376351
    377 !!!******!!
    378 !!! TIME !!
    379 !!!******!!
    380 IF (JULYR .ne. 9999) THEN
    381   !
    382   ! specified
    383   !
    384   ptime = (GMT + elaps/3700.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT
    385   ptime = MODULO(ptime,24.)   !! the two arguments of MODULO must be of the same type
    386   ptime = ptime / 24.
    387   pday = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
    388   pday = MODULO(int(pday),669)
    389   MY = (JULYR-2000) + (88800.*(JULDAY - 1)+3700.*GMT+elaps)/59496000.
    390   MY = INT(MY)
    391 ELSE
    392   !
    393   ! idealized
    394   !
     352!!!***********!!
     353!!! IDEALIZED !!
     354!!!***********!!
     355IF (JULYR .eq. 9999) THEN
    395356  PRINT *,'** Mars ** IDEALIZED SIMULATION'
    396357  PRINT *,'** Mars ** BEWARE: input_coord must be here'
    … …  
    402363  read(14,*) lct_input
    403364  close(14)
    404   ptime = lct_input - lon_input / 15. + elaps/3700.
    405   ptime = MODULO(ptime,24.)
    406   ptime = ptime / 24.
    407   pday = floor(ls2sol(ls_input)) + INT((3700*(lct_input - lon_input / 15.) + elaps)/88800)
    408   pday = MODULO(int(pday),669)
    409   MY = 2024
    410   !day_ini = floor(ls2sol(ls_input)) !! pday at firstcall is day_ini
    411 ENDIF
    412 print *,'** Mars ** TIME IS', pday, ptime*24.
    413 
     365ENDIF
    414366
    415367!----------!
    416368! ALLOCATE !
    417369!----------!
     370!-------------------------------------------------------------------------------!
     371! outputs:                                                                      !       
     372!    pdu(ngrid,nlayermx)        \                                               !
     373!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding    !
     374!    pdt(ngrid,nlayermx)         /  variables due to physical processes.        !
     375!    pdq(ngrid,nlayermx)        /                                               !
     376!    pdpsrf(ngrid)             /                                                !
     377!    tracerdyn                 call tracer in dynamical part of GCM ?           !
     378!-------------------------------------------------------------------------------!
    418379ALLOCATE(pdpsrf(ngrid))
    419380ALLOCATE(pdu(ngrid,nlayer))
    … …  
    450411CALL allocate_comm_wrf(ngrid,nlayer)
    451412#endif
    452 ALLOCATE(q2_val(nlayer+1))
    453 ALLOCATE(qsurf_val(nq))
    454 ALLOCATE(tsoil_val(nsoil))
    455 ALLOCATE(isoil_val(nsoil))
    456 ALLOCATE(dsoil_val(nsoil))
    457413ALLOCATE(pplev(ngrid,nlayer+1))  !!!!!
    458414ALLOCATE(pplay(ngrid,nlayer))    !!!!!
    … …  
    481437!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    482438
    483 !! INITIALIZE AND ALLOCATE EVERYTHING !!
     439#ifndef NOPHYS
     440!! INITIALIZE AND ALLOCATE EVERYTHING !! here, only firstcall
    484441allocation_firstcall: IF (firstcall .EQV. .true.) THEN
    485 
    486   !! TRACERS
    487   ALLOCATE(noms(nq)) !! est fait dans initracer normalement
    488                      !! tableau dans tracer_mod.F90
    489 
    490   !!! name of tracers to mimic entries in tracer.def
    491   !!! ----> IT IS IMPORTANT TO KEEP SAME ORDERING AS IN THE REGISTRY !!!!
    492   !!!                               SAME NAMING   AS IN THE LMD PHYSICS !!!!
    493   SELECT CASE (MARS_MODE)
    494     CASE(0,10) 
    495       noms(nq) = 'co2'
    496     CASE(1)    ! scalar:qh2o,qh2o_ice
    497       noms(1)  = 'h2o_vap'
    498       noms(2)  = 'h2o_ice'
    499     CASE(2)    ! scalar:qdust
    500       noms(1)  = 'dust01'
    501     CASE(3)    ! scalar:qdust,qdustn
    502       noms(1)  = 'dust_mass'
    503       noms(2)  = 'dust_number'
    504     CASE(11)   ! scalar:qh2o,qh2o_ice,qdust,qdustn
    505       noms(1)  = 'h2o_vap'
    506       noms(2)  = 'h2o_ice'
    507       noms(3)  = 'dust_mass'
    508       noms(4)  = 'dust_number'
    509     CASE(12)
    510       noms(1)  = 'h2o_vap'
    511       noms(2)  = 'h2o_ice'
    512       noms(3)  = 'dust_mass'
    513       noms(4)  = 'dust_number'
    514       noms(5)  = 'ccn_mass'
    515       noms(6)  = 'ccn_number'
    516     CASE(20)
    517       noms(1) = 'qtrac1'
    518     CASE(42)
    519       noms(1)  = 'co2'
    520       noms(2)  = 'co'
    521       noms(3)  = 'o'
    522       noms(4)  = 'o1d'
    523       noms(5)  = 'o2'
    524       noms(6)  = 'o3'
    525       noms(7)  = 'h'
    526       noms(8)  = 'h2'
    527       noms(9)  = 'oh'
    528       noms(10)  = 'ho2'
    529       noms(11)  = 'h2o2'
    530       noms(12)  = 'ch4'
    531       noms(13)  = 'n2'
    532       noms(14)  = 'ar'
    533       noms(15)  = 'h2o_ice'
    534       noms(16)  = 'h2o_vap'
    535       noms(17)  = 'dust_mass'
    536       noms(18)  = 'dust_number'
    537   END SELECT
    538 
    539 !!*******************************************!!
    540 
    541 #ifndef NOPHYS
     442  !! tracers' name
     443  PRINT *,'** Mars ** TRACERS NAMES'
     444  CALL update_inputs_physiq_tracers(nq,MARS_MODE)
    542445  !! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM)
    543446  !! parameters are defined in the module_model_constants.F WRF routine
    … …  
    545448  CALL phys_state_var_init(ngrid,nlayer,nq,&
    546449                           wdaysec,ptimestep,1/reradius,g,r_d,cp)
    547 
    548450  !! Fill planetary parameters in modules
    549451  !! Values defined in the module_model_constants.F WRF routine
    550 
    551   !! comcstfi_h
    552   omeg = womeg               
    553   mugaz = wmugaz 
    554   print*,"check: omeg,mugaz"
    555   print*,omeg,mugaz
    556   !! planet_h
    557   year_day = wyear_day
    558   periheli = wperiheli
    559   aphelie = waphelie
    560   peri_day = wperi_day
    561   obliquit = wobliquit
    562   emin_turb = wemin_turb
    563   lmixmin = wlmixmin
    564   print*,"check: year_day,periheli,aphelie,peri_day,obliquit"
    565   print*,year_day,periheli,aphelie,peri_day,obliquit
    566   print*,"check: emin_turb,lmixmin"
    567   print*,emin_turb,lmixmin
    568   !! surfdat_h
    569   emissiv = wemissiv
    570   emisice(1) = wemissiceN
    571   emisice(2) = wemissiceS
    572   albedice(1) = walbediceN
    573   albedice(2) = walbediceS
    574   iceradius(1) = wiceradiusN
    575   iceradius(2) = wiceradiusS
    576   dtemisice(1) = wdtemisiceN
    577   dtemisice(2) = wdtemisiceS
    578   z0_default = wz0
    579   print*,"check: z0def,emissiv,emisice,albedice,iceradius,dtemisice"
    580   print*,z0_default,emissiv,emisice,albedice,iceradius,dtemisice
    581   !! comsoil_h
    582   volcapa = wvolcapa
    583   print*,"check: volcapa",volcapa
    584 
     452  CALL update_inputs_physiq_constants
    585453  !! Read callphys.def !!
    586454  call conf_phys(ngrid,nlayer,nq)
    587 
    588 #endif
    589 
    590455ENDIF allocation_firstcall
     456#endif
    591457
    592458!!*****************************!!
    … …  
    683549ENDDO
    684550
    685 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    686 !!! ONE DOMAIN CASE
    687 !!! --> we simply need to initialize static and physics inputs
    688 !!! SEVERAL DOMAINS
    689 !!! --> we update static and physics inputs each time
    690 !!!     to account for
    691 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    692 pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN
    693 
    694 PRINT *,'** Mars ** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.'
    695 
    696 !!!!!!!!!!!!
    697 !!! GEOM !!!
    698 !!!!!!!!!!!!
    699 
    700 ALLOCATE(plon(ngrid))
    701 ALLOCATE(plat(ngrid))
    702 ALLOCATE(parea(ngrid))
    703 
    704 DO j = jps,jpe
    705 DO i = ips,ipe
    706 
    707   !-----------------------------------!
    708   ! 1D subscript for physics "cursor" !
    709   !-----------------------------------!
    710   subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
    711 
    712   !----------------------------------------!
    713   ! Surface of each part of the grid (m^2) !
    714   !----------------------------------------!
    715   !parea(subs) = dx*dy                           !! 1. idealized cases - computational grid
    716   parea(subs) = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
    717   !parea(subs)=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
    718 
    719   !---------------------------------------------!
    720   ! Mass-point latitude and longitude (radians) !
    721   !---------------------------------------------!
    722   IF (JULYR .ne. 9999) THEN
    723    plat(subs) = XLAT(i,j)*DEGRAD
    724    plon(subs) = XLONG(i,j)*DEGRAD
    725   ELSE
    726    !!! IDEALIZED CASE
    727    IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lat: ',lat_input
    728    IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lon: ',lon_input
    729    plat(subs) = lat_input*DEGRAD
    730    plon(subs) = lon_input*DEGRAD
    731   ENDIF
    732 
    733 ENDDO
    734 ENDDO
    735 
    736 !! FILL GEOMETRICAL ARRAYS !!
    737 call ini_fillgeom(ngrid,plat,plon,parea)
    738 
    739 DEALLOCATE(plon)
    740 DEALLOCATE(plat)
    741 DEALLOCATE(parea)
    742 
    743 !!*******************************************!!
    744 !!*******************************************!!
    745 
    746 DO j = jps,jpe
    747 DO i = ips,ipe
    748 
    749 !-----------------------------------!
    750 ! 1D subscript for physics "cursor" !
    751 !-----------------------------------!
    752 subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
    753 
    754 !! **************
    755 !! use slope_mod
    756 !! **************
    757 !-------------------!
    758 ! Slope inclination !
    759 !-------------------!
    760 IF (JULYR .ne. 9999) THEN
    761   theta_val=atan(sqrt( (1000.*SLPX(i,j))**2 + (1000.*SLPY(i,j))**2 ))
    762   theta_val=theta_val/DEGRAD
    763 ELSE
    764   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope inclination is 0'
    765   theta_val=0.
    766 ENDIF
    767 !-------------------------------------------!
    768 ! Slope orientation; 0 is north, 90 is east !
    769 !-------------------------------------------!
    770 IF (JULYR .ne. 9999) THEN
    771   psi_val=-90.*DEGRAD-atan(SLPY(i,j)/SLPX(i,j))
    772   if (SLPX(i,j) .ge. 0.) then
    773     psi_val=psi_val-180.*DEGRAD
    774   endif
    775   psi_val=360.*DEGRAD+psi_val
    776   psi_val=psi_val/DEGRAD
    777   psi_val = MODULO(psi_val+180.,360.)
    778 ELSE
    779   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope orientation is 0 (well, whatever)'
    780   psi_val=0.
    781 ENDIF
    782 theta_sl(subs) = theta_val
    783 psi_sl(subs) = psi_val
    784 
    785 !! **************
    786 !! use surfdat_h
    787 !! **************
    788551!---------------------------------------------------------!
    789552! Ground geopotential                                     !
    … …  
    791554!---------------------------------------------------------!
    792555phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.)   !! NB: z_prof are half levels, so z_prof(1) is not the surface
    793 !---------------------------------!
    794 ! Ground albedo & Thermal Inertia !
    795 !---------------------------------!
    796 IF (JULYR .ne. 9999) THEN
    797  IF (CST_AL == 0) THEN
    798  albedodat_val=M_ALBEDO(i,j)
    799  ELSE
    800  albedodat_val=CST_AL
    801  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT ALBEDO ', albedodat_val
    802  ENDIF
    803 ELSE
    804  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION albedo: ', CST_AL
    805  albedodat_val=CST_AL
    806 ENDIF
    807 !-----------------------------------------!
    808 ! Gravity wave parametrization            !
    809 ! NB: usually 0 in mesoscale applications !
    810 !-----------------------------------------!
    811 IF (JULYR .ne. 9999) THEN
    812  zmea_val=M_GW(i,1,j)
    813  zstd_val=M_GW(i,2,j)
    814  zsig_val=M_GW(i,3,j)
    815  zgam_val=M_GW(i,4,j)
    816  zthe_val=M_GW(i,5,j)
    817 ELSE
    818  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION GWdrag OFF'
    819  zmea_val=0.
    820  zstd_val=0.
    821  zsig_val=0.
    822  zgam_val=0.
    823  zthe_val=0.
    824 ENDIF
    825 phisfi(subs) = phisfi_val
    826 albedodat(subs) = albedodat_val
    827 zmea(subs) = zmea_val
    828 zstd(subs) = zstd_val
    829 zsig(subs) = zsig_val
    830 zgam(subs) = zgam_val
    831 zthe(subs) = zthe_val
    832 !----------------------------!
    833 ! Variable surface roughness !
    834 !----------------------------!
    835 IF (JULYR .ne. 9999) THEN
    836  IF (CST_Z0 == 0) THEN
    837    z0_val = M_Z0(i,j)
    838  ELSE
    839    z0_val = CST_Z0
    840    IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
    841  ENDIF
    842 ELSE
    843  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
    844  z0_val=CST_Z0
    845 ENDIF
    846 !!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
    847 IF (z0_val == 0.) THEN
    848    IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
    849    z0_val = 0.01
    850 ENDIF
    851 !!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
    852 IF (z0_val < 0.) THEN
    853    PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
    854    PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
    855    PRINT *, '           -- or check the constant value set in namelist.input'
    856    STOP
    857 ENDIF
    858 z0(subs) = z0_val
    859 !-----------------------------------------------!
    860 ! Ground temperature, emissivity, CO2 ice cover !
    861 !-----------------------------------------------!
    862 tsurf(subs) = M_TSURF(i,j)
    863 emis(subs) = M_EMISS(i,j)
    864 co2ice(subs) = M_CO2ICE(i,j)
    865 !-------------------!
    866 ! Tracer at surface !
    867 !-------------------!
    868 qsurf_val(:)=0. ! default case
    869 SELECT CASE (MARS_MODE)
    870     CASE(1)
    871     qsurf_val(2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
    872                                    !! ----- retrocompatible ancienne physique
    873                                    !! ----- [H2O ice is last tracer in qsurf in LMD physics]
    874     CASE(2)   
    875     qsurf_val(1)=0.                !! not coupled with lifting for the moment [non remobilise]
    876     CASE(3)
    877     qsurf_val(1)=q_prof(1,1)                !!! temporaire, a definir       
    878     qsurf_val(2)=q_prof(1,2)     
    879     CASE(11)
    880     qsurf_val(2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
    881     qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
    882     CASE(12)
    883     qsurf_val(2)=M_H2OICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
    884     qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
    885 END SELECT
    886 qsurf(subs,:) = qsurf_val(:)
    887 
    888 !! **************
    889 !! use comsoil_h
    890 !! **************
    891 !---------------------------------!
    892 ! Ground albedo & Thermal Inertia !
    893 !---------------------------------!
    894 IF (JULYR .ne. 9999) THEN
    895  IF (CST_TI == 0) THEN
    896  inertiedat_val=M_TI(i,j)
    897  ELSE
    898  inertiedat_val=CST_TI
    899  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT THERMAL INERTIA ', inertiedat_val
    900  ENDIF
    901 ELSE
    902  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION inertia: ',CST_TI
    903  inertiedat_val=CST_TI
    904 ENDIF
    905 !inertiedat(subs) = inertiedat_val
    906 !--pb de dimensions???!!???
    907 IF (JULYR .ne. 9999) THEN
    908   isoil_val(:)=M_ISOIL(i,:,j)
    909   dsoil_val(:)=M_DSOIL(i,:,j)
    910 ELSE
    911    IF ( nsoil .lt. 18 ) THEN
    912        PRINT *,'** Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
    913        STOP
    914    ENDIF
    915    IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
    916   do k=1,nsoil
    917    isoil_val(k) = inertiedat_val
    918 !  dsoil_val(k) = sqrt(887.75/3.14)*((2.**(k-0.5))-1.) * inertiedat_val / wvolcapa    ! old setting
    919    dsoil_val(k) = 2.E-4 * (2.**(k-0.5-1.))                                            ! new gcm settings
    920    IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** ISOIL DSOIL are ',isoil_val(k), dsoil_val(k)
    921   enddo
    922 ENDIF
    923 inertiedat(subs,:) = isoil_val(:)
    924 mlayer(0:nsoil-1)=dsoil_val(:)
    925     !!!!!!!!!!!!!!!!! DONE in soil_setting.F
    926     ! 1.5 Build layer(); following the same law as mlayer()
    927     ! Assuming layer distribution follows mid-layer law:
    928     ! layer(k)=lay1*alpha**(k-1)
    929     lay1=sqrt(mlayer(0)*mlayer(1))
    930     alpha=mlayer(1)/mlayer(0)
    931     do iii=1,nsoil
    932        layer(iii)=lay1*(alpha**(iii-1))
    933     enddo
    934 !------------------------!
    935 ! Deep soil temperatures !
    936 !------------------------!
    937 IF (M_TSOIL(i,1,j) .gt. 0. .and. JULYR .ne. 9999) THEN
    938   tsoil_val(:)=M_TSOIL(i,:,j)
    939 ELSE
    940   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** no tsoil. set it to tsurf.'
    941   do k=1,nsoil
    942    tsoil_val(k) = M_TSURF(i,j)
    943   enddo
    944 ENDIF
    945 tsoil(subs,:) = tsoil_val(:)
    946 
    947 !! **************
    948 !! use turb_mod
    949 !! **************
    950 IF (.not. restart) THEN
    951    q2_val(:) = 1.E-6      !PBL wind variance
    952    wstar_val=0.
    953 ELSE
    954    q2_val(:)=M_Q2(i,:,j)
    955    wstar_val=M_WSTAR(i,j)
    956 ENDIF
    957 q2(subs,:) = q2_val(:)
    958 wstar(subs) = wstar_val
    959 
    960 !! ********************
    961 !! use dimradmars_mod
    962 !! ********************
    963 IF (.not. restart) THEN
    964    fluxrad_val=0.
    965 ELSE
    966    fluxrad_val=M_FLUXRAD(i,j)
    967 ENDIF
    968 fluxrad(subs) = fluxrad_val
    969 !! et fluxrad_sky ???!???
    970 
    971 ENDDO
    972 ENDDO
    973 
    974 !!---------------------!!
    975 !! OUTPUT FOR CHECKING !!
    976 !!---------------------!!
    977 print*,"check: theta_sl",theta_sl(1),theta_sl(ngrid)
    978 print*,"check: psi_sl",psi_sl(1),psi_sl(ngrid)
    979 print*,"check: phisfi",phisfi(1),phisfi(ngrid)
    980 print*,"check: albedodat",albedodat(1),albedodat(ngrid)
    981 print*,"check: zmea",zmea(1),zmea(ngrid)
    982 print*,"check: zstd",zstd(1),zstd(ngrid)
    983 print*,"check: zsig",zsig(1),zsig(ngrid)
    984 print*,"check: zgam",zgam(1),zgam(ngrid)
    985 print*,"check: zthe",zthe(1),zthe(ngrid)
    986 print*,"check: z0",z0(1),z0(ngrid)
    987 print*,"check: tsurf",tsurf(1),tsurf(ngrid)
    988 print*,"check: emis",emis(1),emis(ngrid)
    989 print*,"check: co2ice",co2ice(1),co2ice(ngrid)
    990 print*,"check: qsurf",qsurf(1,:),qsurf(ngrid,:)
    991 print*,"check: inertiedat",inertiedat(1,:),inertiedat(ngrid,:)
    992 print*,"check: mlayer",mlayer(:)
    993 print*,"check: layer",layer(:)
    994 print*,"check: q2",q2(1,:),q2(ngrid,:)
    995 print*,"check: wstar",wstar(1),wstar(ngrid)
    996 print*,"check: fluxrad",fluxrad(1),fluxrad(ngrid)
    997 print*,"check: tsoil",tsoil(1,:),tsoil(ngrid,:)
    998 print*,"check: noms",noms
    999 
    1000 ENDIF pass_interface
    1001556
    1002557!!*****************!!
    1003558!! CLEAN THE PLACE !!
    1004559!!*****************!!
    1005 DEALLOCATE(q2_val)
    1006 DEALLOCATE(qsurf_val)
    1007 DEALLOCATE(tsoil_val)
    1008 DEALLOCATE(isoil_val)
    1009 DEALLOCATE(dsoil_val)
    1010560DEALLOCATE(dz8w_prof)
    1011561DEALLOCATE(z_prof)
    … …  
    1018568DEALLOCATE(q_prof)
    1019569
     570
     571!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     572!!! ONE DOMAIN CASE
     573!!! --> we simply need to initialize static and physics inputs
     574!!! SEVERAL DOMAINS
     575!!! --> we update static and physics inputs each time
     576!!!     to account for domain change
     577!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     578pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN
     579PRINT *,'** Mars ** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.'
     580!!*******************************************!!
     581!!*******************************************!!
     582CALL update_inputs_physiq_geom( &
     583            ims,ime,jms,jme,&
     584            ips,ipe,jps,jpe,&
     585            JULYR,ngrid,&
     586            DX,DY,MSFT,&
     587            lat_input, lon_input,&
     588            XLAT,XLONG)
     589!!!
     590CALL update_inputs_physiq_slope( &
     591            ims,ime,jms,jme,&
     592            ips,ipe,jps,jpe,&
     593            JULYR,&
     594            SLPX,SLPY)
     595!!!
     596CALL update_inputs_physiq_soil( &
     597            ims,ime,jms,jme,&
     598            ips,ipe,jps,jpe,&
     599            JULYR,nsoil,&
     600            M_TI,CST_TI,&
     601            M_ISOIL,M_DSOIL,&
     602            M_TSOIL,M_TSURF)
     603!!!
     604CALL update_inputs_physiq_surf( &
     605            ims,ime,jms,jme,&
     606            ips,ipe,jps,jpe,&
     607            JULYR,MARS_MODE,&
     608            M_ALBEDO,CST_AL,&
     609            M_TSURF,M_EMISS,M_CO2ICE,&
     610            M_GW,M_Z0,&
     611            M_H2OICE,&
     612            phisfi_val)
     613!!!
     614CALL update_inputs_physiq_turb( &
     615            ims,ime,jms,jme,kms,kme,&
     616            ips,ipe,jps,jpe,&
     617            RESTART,isles,&
     618            M_Q2,M_WSTAR)
     619!!!
     620CALL update_inputs_physiq_rad( &
     621            ims,ime,jms,jme,&
     622            ips,ipe,jps,jpe,&
     623            RESTART,&
     624            M_FLUXRAD)
     625!!!
     626ENDIF pass_interface
     627!!*******************************************!!
     628!!*******************************************!!
     629
    1020630!!!!!!!!!!!!!!!!!!!!!!
    1021631!!!!!!!!!!!!!!!!!!!!!!
    … …  
    1025635
    1026636call_physics : IF (wappel_phys .ne. 0.) THEN
    1027 
    1028 !-------------------------------------------------------------------------------!
    1029 ! outputs:                                                                      !       
    1030 !    pdu(ngrid,nlayermx)        \                                               !
    1031 !    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding    !
    1032 !    pdt(ngrid,nlayermx)         /  variables due to physical processes.        !
    1033 !    pdq(ngrid,nlayermx)        /                                               !
    1034 !    pdpsrf(ngrid)             /                                                !
    1035 !    tracerdyn                 call tracer in dynamical part of GCM ?           !
    1036 !-------------------------------------------------------------------------------!
     637!!! initialize tendencies (security)
    1037638pdpsrf(:)=0.
    1038639pdu(:,:)=0.
    … …  
    1042643#ifndef NOPHYS
    1043644print *, '** Mars ** CALL TO LMD PHYSICS'
     645!!!
     646CALL update_inputs_physiq_time(&
     647            JULYR,JULDAY,GMT,&
     648            elaps,&
     649            lct_input,lon_input,ls_input,&
     650            ptime,pday,MY)
     651!!!
    1044652CALL physiq (ngrid,nlayer,nq,                         &
    1045653             firstcall,lastcall,pday,ptime,ptimestep, &
    … …  
    1083691
    1084692!! OUTPUT OUTPUT OUTPUT
     693!-------------------------------------------------------!
     694! Save key variables for restart and output and nesting ! 
     695!-------------------------------------------------------!
     696!!!
     697CALL update_outputs_physiq_surf( &
     698            ims,ime,jms,jme,&
     699            ips,ipe,jps,jpe,&
     700            MARS_MODE,&
     701            M_TSURF,M_CO2ICE,&
     702            M_H2OICE)
     703!!!
     704CALL update_outputs_physiq_soil( &
     705            ims,ime,jms,jme,&
     706            ips,ipe,jps,jpe,&
     707            nsoil,&
     708            M_TSOIL)
     709!!!
     710CALL update_outputs_physiq_rad( &
     711            ims,ime,jms,jme,&
     712            ips,ipe,jps,jpe,&
     713            M_FLUXRAD)
     714!!!
     715CALL update_outputs_physiq_turb( &
     716            ims,ime,jms,jme,kms,kme,&
     717            ips,ipe,jps,jpe,&
     718            M_Q2,M_WSTAR,&
     719            HFMAX,ZMAX,USTM,HFX)
     720
    1085721DO j = jps,jpe
    1086722DO i = ips,ipe
    1087723
    1088724  subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
    1089 
    1090   !-------------------------------------------------------!
    1091   ! Save key variables for restart and output and nesting ! 
    1092   !-------------------------------------------------------!
    1093   M_TSOIL(i,:,j) = tsoil(subs,:)
    1094   M_CO2ICE(i,j) = co2ice(subs)
    1095   M_Q2(i,kps:kpe+1,j) = q2(subs,:)
    1096   M_TSURF(i,j) = tsurf(subs)
    1097   M_WSTAR(i,j) = wstar(subs)
    1098   M_FLUXRAD(i,j) = fluxrad(subs)
    1099   SELECT CASE (MARS_MODE)
    1100    CASE (1,11,12)
    1101      M_H2OICE(i,j) = qsurf(subs,2)  !! see above Tracer at surface
    1102   END SELECT
    1103 
    1104   !! output only (arrays already in phys modules)
    1105   HFMAX(i,j) = HFMAX_TH(subs)
    1106   ZMAX(i,j) = ZMAX_TH(subs)
    1107   USTM(i,j) = ustar(subs)
    1108   HFX(i,j) = sensibFlux(subs)
    1109725
    1110726  !! output only (cf comm_wrf)
    … …  
    1200816END SUBROUTINE lmd_driver
    1201817
    1202 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1203       real function ls2sol(ls)
    1204 
    1205 !c  Returns solar longitude, Ls (in deg.), from day number (in sol),
    1206 !c  where sol=0=Ls=0 at the northern hemisphere spring equinox
    1207 
    1208       implicit none
    1209 
    1210 !c  Arguments:
    1211       real ls
    1212 
    1213 !c  Local:
    1214       double precision xref,zx0,zteta,zz
    1215 !c      xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
    1216       double precision year_day
    1217       double precision peri_day,timeperi,e_elips
    1218       double precision pi,degrad
    1219       parameter (year_day=668.6d0) ! number of sols in a amartian year
    1220 !c      data peri_day /485.0/
    1221       parameter (peri_day=485.35d0) ! date (in sols) of perihelion
    1222 !c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
    1223       parameter (timeperi=1.90258341759902d0)
    1224       parameter (e_elips=0.0934d0)  ! eccentricity of orbit
    1225       parameter (pi=3.14159265358979d0)
    1226       parameter (degrad=57.2957795130823d0)
    1227 
    1228       if (abs(ls).lt.1.0e-5) then
    1229          if (ls.ge.0.0) then
    1230             ls2sol = 0.0
    1231          else
    1232             ls2sol = year_day
    1233          end if
    1234          return
    1235       end if
    1236 
    1237       zteta = ls/degrad + timeperi
    1238       zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
    1239       xref = zx0-e_elips*dsin(zx0)
    1240       zz = xref/(2.*pi)
    1241       ls2sol = zz*year_day + peri_day
    1242       if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day
    1243       if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day
    1244 
    1245       return
    1246       end function ls2sol
    1247 !!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1248 
    1249818END MODULE module_lmd_driver
Note: See TracChangeset for help on using the changeset viewer.