Changeset 2868


Ignore:
Timestamp:
Jan 16, 2023, 4:47:08 PM (23 months ago)
Author:
jleconte
Message:

Big cleanup of interface_v4. M_TSURF is now P_TSURF.

Location:
trunk/WRF.COMMON/INTERFACES_V4
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/callphysiq_mod.F

    r2866 r2868  
    4343
    4444  call allocate_comm_wrf(klon,llm)
    45 
    4645! Call physics package with required inputs/outputs
    4746  CALL physiq(klon,           & ! ngrid
    4847              llm,            & ! nlayer
    4948              nqtot,          & ! nq
    50 !              noms,          & ! nametrac
    5149              debut_split,    & ! firstcall
    5250              lafin_split,    & ! lastcall

  • trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_inputs_physiq_mod.F

    r2866 r2868  
    2727  REAL :: sec,nsec
    2828 
    29   IF (JULYR .le. 8999) THEN
    30     if (tlocked .eqv. .false.) THEN
     29  !IF (JULYR .le. 8999) THEN
     30  !  if (tlocked .eqv. .false.) THEN
    3131      JH_cur_split = (GMT + elaps/3600.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
    3232      JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type
    … …  
    3636      MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600.0*GMT+elaps)/31968000
    3737      MY = INT(MY)
    38     ELSE
    39       JH_cur_split = (GMT)! + elaps/420000.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
    40       JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type
    41       JH_cur_split = JH_cur_split / 24.
    42       JD_cur = (JULDAY - 1 + INT((3600*GMT)/86400))!+elaps)/1.008e7))
    43       JD_cur = MODULO(int(JD_cur),2)
    44       MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600*GMT+elaps)/31968000
    45       MY = INT(MY)
    46     ENDIF
    47   ELSE
    48     if (tlocked .eqv. .false.) THEN
    49       JH_cur_split = lct_input - lon_input / 15. + elaps/1500.0
    50       JD_cur =  INT((sec*(lct_input - lon_input / 15.) + elaps)/36000)
    51       !ptime = lct_input - lon_input / 15. + elaps/3600.
    52       !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
    53     ELSE
    54       JH_cur_split = lct_input - lon_input / 15. !+ elaps/1500.0
    55       !pday =  INT((sec*(lct_input - lon_input / 15.)+ elaps)/36000)
    56       JD_cur =  INT((sec*(lct_input - lon_input / 15.))/3600)
    57       !print*,'ptime',ptime
    58       !print*,'pday',pday
    59     !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
    60       JH_cur_split = MODULO(ptime,24.)
    61       JH_cur_split = JH_cur_split / 24.
    62       JD_cur = MODULO(int(pday),365)
    63       MY = 2024
    64     ENDIF
    65   ENDIF
     38  !  ELSE
     39  !    JH_cur_split = (GMT)! + elaps/420000.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT
     40  !    JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type
     41  !    JH_cur_split = JH_cur_split / 24.
     42  !    JD_cur = (JULDAY - 1 + INT((3600*GMT)/86400))!+elaps)/1.008e7))
     43  !    JD_cur = MODULO(int(JD_cur),2)
     44  !    MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600*GMT+elaps)/31968000
     45  !    MY = INT(MY)
     46  !  ENDIF
     47  !ELSE
     48  !  if (tlocked .eqv. .false.) THEN
     49  !    JH_cur_split = lct_input - lon_input / 15. + elaps/1500.0
     50  !    JD_cur =  INT((sec*(lct_input - lon_input / 15.) + elaps)/36000)
     51  !    !ptime = lct_input - lon_input / 15. + elaps/3600.
     52  !    !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
     53  !  ELSE
     54  !    JH_cur_split = lct_input - lon_input / 15. !+ elaps/1500.0
     55  !    !pday =  INT((sec*(lct_input - lon_input / 15.)+ elaps)/36000)
     56  !    JD_cur =  INT((sec*(lct_input - lon_input / 15.))/3600)
     57  !    !print*,'ptime',ptime
     58  !    !print*,'pday',pday
     59  !  !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
     60  !    JH_cur_split = MODULO(ptime,24.)
     61  !    JH_cur_split = JH_cur_split / 24.
     62  !    JD_cur = MODULO(int(pday),365)
     63  !    MY = 2024
     64  !  ENDIF
     65  !ENDIF
     66
    6667
    6768END SUBROUTINE update_inputs_physiq_time
    … …  
    8081     
    8182                 !! tableau dans tracer_mod.F90
    82   nqtot=nq
    8383  if ((TRACER_MODE .eq. 1) .or. (TRACER_MODE .ge. 42)) THEN
    8484    nqtot=2
    … …  
    107107SUBROUTINE update_inputs_physiq_constants
    108108
    109    !USE module_model_constants
    110    !use comcstfi_h, only: omeg,mugaz
    111    !use planete_h,  only: year_day,periheli,aphelie, &
    112    !                       peri_day,obliquit,emin_turb, &
    113    !                       lmixmin
    114109   use planete_mod, only: year_day, periastr, apoastr, peri_day,&
    115110                       obliquit, z0, lmixmin, emin_turb
    … …  
    117112                         emisice,dtemisice
    118113   !                       z0_default
    119    !use comsoil_h, only: volcapa
    120114   use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
    121    !! comcstfi_h
    122115   use phys_state_var_mod, only :cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tsea_ice
    123    !use iniorbit
    124    !use time_phylmdz_mod, only: dtphys, daysec,day_ini
    125  
    126  
    127      !open(17,file='controle.txt',form='formatted',status='old')
    128      !rewind(17)
    129      !read(17,*)
    130      !read(17,*)
    131      !read(17,*) day_ini !(tab0+3)
    132      !read(17,*)
    133      !read(17,*) !tab0+5)
    134      !read(17,*) omeg !(tab0+6)
    135      !read(17,*) !(tab0+7)
    136      !read(17,*)  !(tab0+8)
    137      !read(17,*)  !(tab0+9)
    138      !read(17,*) daysec
    139      !read(17,*) dtphys !tab0+11)
    140      !read(17,*)
    141      !read(17,*)
    142      !read(17,*) year_day !(tab0+14)
    143      !read(17,*) periastr !tab0+15)
    144      !read(17,*) apoastr !tab0+16)
    145      !read(17,*) peri_day !tab0+17)
    146      !read(17,*) obliquit !tab0+18)
    147      !read(17,*) z0
    148      !read(17,*)
    149      !read(17,*)
    150      !read(17,*)
    151      !read(17,*)
    152      !read(17,*) emisice(1)
    153      !read(17,*) emisice(2)
    154      !read(17,*) emissiv
    155      !read(17,*)
    156      !read(17,*)
    157      !read(17,*)
    158      !read(17,*)
    159      !read(17,*) iceradius(1)
    160      !read(17,*) iceradius(2)
    161      !read(17,*) dtemisice(1)
    162      !read(17,*) dtemisice(2)
    163      !close(17)
    164      !cpp=(8.314511/(mugaz/1000.0))/rcp
    165      !print*,'cpp',cpp
    166      !print*,'g',g
    167 
    168      !emissiv(:)=EMIS
    169      !cloudfrac(:,:)=0.5
    170      !totcloudfrac(:)=0.5
    171      !hice(:)=0.
    172      !rnat(:)=0.
    173      !pctsrf_sic(:)=0.
    174      !tsea_ice(:)=0.
    175      !qsurf(:,:) = 0.
    176      !print*,'iceradius',iceradius,'dtemisice',dtemisice
    177      !print*,'apoastr,periastr,year_day,peri_day,obliq',apoastr,periastr,year_day,peri_day,obliquit
    178      !print*,'emissiv',emissiv
    179  
    180      ! SUBROUTINE iniorbit(apoastr,periastr,year_day,peri_day,obliq)
    181 
     116   !JL22 this routine does not do anything for the generic interface
     117   ! The various use abave can surely be removed.
    182118END SUBROUTINE update_inputs_physiq_constants
    183119
    … …  
    272208   latitude_deg(:) = plat(:)/DEGRAD
    273209   cell_area(:) = parea(:)
    274    !call planetwide_sumval(parea,totarea_planet)
    275    !print*,'parea',parea(1)
    276    !totarea=SSUM(ngrid,parea,1)
    277210   totarea=ngrid*parea(1)
    278    !totarea_planet=SSUM(ngrid,parea,1)
    279211   totarea_planet=ngrid*parea(1)
    280212
    … …  
    294226   DO k=1, nlayer
    295227   read(12,*) znw(k)
    296    !write(6,*) 'read level ', k,grid%znw(k)
    297228   ENDDO
    298229   close(12)
    … …  
    319250            JULYR,TRACER_MODE,&
    320251            M_ALBEDO,CST_AL,&
    321             M_TSURF,M_EMISS,M_CO2ICE,&
     252            P_TSURF,M_EMISS,M_CO2ICE,&
    322253            M_GW,M_Z0,CST_Z0,&
    323254            M_H2OICE,&
    … …  
    335266   REAL, INTENT(IN  ) :: CST_AL, phisfi_val, CST_Z0
    336267   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
    337      M_ALBEDO,M_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
     268     M_ALBEDO,P_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
    338269   REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: M_GW 
    339270
    340    !print*,'ALLOCATED(phisfi)',ALLOCATED(phisfi)
    341    !print*,'size phisfi',size(phisfi)
    342271   DO j = jps,jpe
    343272   DO i = ips,ipe
    … …  
    352281     !---------------------!
    353282     phisfi(subs) = phisfi_val
    354 !     print*,'size phisfi',size(phisfi)
    355      !print*,'phisfi',phisfi(subs)
    356283     !---------------!
    357284     ! Ground albedo !
    … …  
    383310     !----------------------------!
    384311     z0 = CST_Z0
    385      !IF (JULYR .le. 8999) THEN
    386      ! IF (CST_Z0 == 0) THEN
    387      !  z0(subs) = M_Z0(i,j)
    388      ! ELSE
    389      !  z0(subs) = CST_Z0
    390      !  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
    391      ! ENDIF
    392      !ELSE
    393      ! IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
    394      ! z0(subs)=CST_Z0
    395      !ENDIF
    396      !!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
    397      !IF (z0(subs) == 0.) THEN
    398      !  IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
    399      !  z0(subs) = 0.01
    400      !ENDIF
    401      !!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
    402      !IF (z0(subs) < 0.) THEN
    403      !  PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
    404      !  PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
    405      !  PRINT *, '           -- or check the constant value set in namelist.input'
    406      !  STOP
    407      !ENDIF
    408 
     312     
    409313     !-----------------------------------------------!
    410314     ! Ground temperature, emissivity, CO2 ice cover !
    411315     !-----------------------------------------------!
    412      tsurf(subs) = M_TSURF(i,j)
     316     tsurf(subs) = P_TSURF(i,j)
    413317     emis(subs) = M_EMISS(i,j)     
    414318     !do i=1,noceanmx
    … …  
    466370            M_TI,CST_TI,&
    467371            M_ISOIL,M_DSOIL,&
    468             M_TSOIL,M_TSURF)
     372            M_TSOIL,P_TSURF)
    469373
    470374   use comsoil_h, only: inertiedat,mlayer,layer,volcapa
    … …  
    476380   REAL, INTENT(IN  ) :: CST_TI
    477381   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
    478      M_TI, M_TSURF
     382     M_TI, P_TSURF
    479383   REAL, DIMENSION( ims:ime, nsoil, jms:jme ), INTENT(IN)  :: &
    480384     M_TSOIL, M_ISOIL, M_DSOIL
    … …  
    541445       IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** no tsoil. set it to tsurf.'
    542446       do k=1,nsoil
    543         !print*,'M_TSURF(i,j)',M_TSURF(1,:)
    544         !print*,'size M_TSURF',size(M_TSURF)
    545         !print*,'size tsoil',size(tsoil)
    546         tsoil(subs,k) = M_TSURF(i,j)
    547         !print*,'tsoil(subs,k)',tsoil(subs,k)
     447        tsoil(subs,k) = P_TSURF(i,j)
    548448       enddo
    549449     ENDIF

  • trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_outputs_physiq_mod.F

    r2866 r2868  
    99            ips,ipe,jps,jpe,&
    1010            TRACER_MODE,&
    11             M_TSURF,M_CO2ICE,&
     11            P_TSURF,M_CO2ICE,&
    1212            M_H2OICE)
    1313
    … …  
    2020   INTEGER :: i,j,subs
    2121   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
    22      M_TSURF,M_CO2ICE,M_H2OICE
    23 
    24    DO j = jps,jpe
    25    DO i = ips,ipe
    26 
    27      !-----------------------------------!
    28      ! 1D subscript for physics "cursor" !
    29      !-----------------------------------!
    30      subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
    31 
    32      !-------------------------------------------------------!
    33      ! Save key variables for restart and output and nesting ! 
    34      !-------------------------------------------------------!
    35      M_TSURF(i,j) = tsurf(subs)
     22     P_TSURF,M_CO2ICE,M_H2OICE
     23
     24   DO j = jps,jpe
     25   DO i = ips,ipe
     26
     27     !-----------------------------------!
     28     ! 1D subscript for physics "cursor" !
     29     !-----------------------------------!
     30     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
     31
     32     !-------------------------------------------------------!
     33     ! Save key variables for restart and output and nesting ! 
     34     !-------------------------------------------------------!
     35     P_TSURF(i,j) = tsurf(subs)
    3636
    3737   ENDDO

  • trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F

    r2866 r2868  
    2626        P8W,DZ8W,T8W,Z,HT,MUT,DNW, &
    2727        U,V,TH,T,P,EXNER,RHO, &
     28        P_HYD, P_HYD_W, &
    2829        PTOP, &
    2930        RADT, &
    … …  
    3536        planet_type, &
    3637        M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
    37         M_H2OICE,M_TSOIL,M_Q2,M_TSURF, &
     38        M_H2OICE,M_TSOIL,M_Q2,P_TSURF, &
    3839        M_FLUXRAD,M_WSTAR,M_ISOIL,M_DSOIL,&
    3940        M_Z0, CST_Z0, M_GW, &
    … …  
    4142        CST_AL, CST_TI, &
    4243        isfflx, diff_opt, km_opt, &
    43         HISTORY_INTERVAL, &
    44         HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,&
     44        HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,&
    4545        CLOUDFRAC,TOTCLOUDFRAC,RH, &
    4646        DQICE,DQVAP,REEVAP,SURFRAIN,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
    … …  
    9494REAL, INTENT(IN  ) :: CST_AL, CST_TI
    9595REAL, INTENT(IN  ) :: PTOP
    96 INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
    9796! 2D arrays         
    9897REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
    … …  
    103102     SLPX,SLPY, &
    104103     M_CO2ICE,M_H2OICE, &
    105      M_TSURF, M_Z0, &
     104     P_TSURF, M_Z0, &
    106105     M_FLUXRAD,M_WSTAR, &
    107106     PSFC,TSK
    … …  
    116115! 3D arrays
    117116REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
    118      dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th
     117     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th,p_hyd,p_hyd_w
    119118REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
    120119     RTHBLTEN,RUBLTEN,RVBLTEN, &
    121      HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,RDUST,VMR_ICE,RICE,&
     120     HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,RDUST,VMR_ICE,RICE,&
    122121     CLOUDFRAC,RH,DQICE,DQVAP,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF
    123122REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
    … …  
    156155   ! <------ outputs:
    157156   !     physical tendencies
    158    REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdtheta
    159157   ! ... intermediate arrays
    160158   REAL, DIMENSION(:), ALLOCATABLE  :: &
    … …  
    171169   LOGICAL, SAVE :: firstcall = .true.
    172170   INTEGER, SAVE :: previous_id = 0
    173 !**************************************************
    174 ! IMPORTANT: pour travailler avec grid nesting
    175 ! --- deux comportements distincts du save
    176 ! ... ex1: ferme planeto avec PGF+MPI: standard
    177 ! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE
    178 !**************************************************
    179 #ifdef SPECIAL_NEST_SAVE
    180       !  une dimension supplementaire liee au nest
    181       REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
    182              dp_save
    183       REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
    184              du_save, dv_save, dt_save,dtheta_save
    185       REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
    186              dq_save     
    187 #else
    188       REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
    189              dp_save
    190       REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
    191              du_save, dv_save, dt_save,dtheta_save
    192       REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
    193              dq_save     
    194 #endif
     171   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: dp_save
     172   REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: du_save, dv_save, dt_save
     173   REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq_save     
    195174
    196175!!!IDEALIZED IDEALIZED
    … …  
    270249! DIMENSIONS FOR THE PHYSICS !
    271250!----------------------------!
    272 !day_ini = JULDAY - 1      !! GCM convention  !! pday at firstcall is day_ini
    273251wappel_phys = RADT
    274252zdt_split = dt*wappel_phys            ! physical timestep (s)
    … …  
    282260lastcall = .false.
    283261! **** needed but hardcoded
    284 
    285           !PRINT *, ips, ipe, jps, jpe
    286           !PRINT *, ngrid
    287262
    288263elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
    … …  
    300275  !firstcall=.false. !! just in case you'd want to get rid of the physics
    301276test=0
    302 #ifdef SPECIAL_NEST_SAVE
    303 PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
    304 IF( .NOT. ALLOCATED( dp_save  ) ) THEN
    305    PRINT *, '**** check **** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
    306    !! here are the arrays that would be useful to save physics tendencies
    307    ALLOCATE(dp_save(ngrid,max_dom))
    308    ALLOCATE(du_save(ngrid,nlayer,max_dom))
    309    ALLOCATE(dv_save(ngrid,nlayer,max_dom))
    310    ALLOCATE(dt_save(ngrid,nlayer,max_dom))
    311    ALLOCATE(dtheta_save(ngrid,nlayer,max_dom))
    312    ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
    313    dp_save(:,:)=0.    !! initialize these arrays ...
    314    du_save(:,:,:)=0.
    315    dv_save(:,:,:)=0.
    316    dt_save(:,:,:)=0.
    317    dtheta_save(:,:,:)=0.
    318    dq_save(:,:,:,:)=0.
    319 ENDIF
    320 IF (id .lt. max_dom) THEN
    321    flag_first_restart=.true.
    322 ELSE
    323    flag_first_restart=.false.
    324 ENDIF
    325 #else
    326277IF( .NOT. ALLOCATED( dp_save  ) ) THEN
    327278ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
    … …  
    329280ALLOCATE(dv_save(ngrid,nlayer))
    330281ALLOCATE(dt_save(ngrid,nlayer))
    331 ALLOCATE(dtheta_save(ngrid,nlayer))
    332282ALLOCATE(dq_save(ngrid,nlayer,nq))
    333283ENDIF
    … …  
    336286dv_save(:,:)=0.
    337287dt_save(:,:)=0.
    338 dtheta_save(:,:)=0.
    339288dq_save(:,:,:)=0.
    340289flag_first_restart=.false.
    341 #endif
    342 ELSE
     290
     291ELSE ! is_first_step
    343292!-------------------------------------------------!
    344293! what is done for the other steps of simulation  !
    … …  
    396345!----------!
    397346CALL allocate_interface(ngrid,nlayer,nq)
    398 ALLOCATE(pdtheta(ngrid,nlayer))
    399347!!!
    400348!!! BIG LOOP : 1. no call for physics, used saved values
    … …  
    402350IF (test.NE.0) THEN
    403351print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test
    404 #ifdef SPECIAL_NEST_SAVE
    405 zdpsrf_omp(:)=dp_save(:,id)
    406 zdufi_omp(:,:)=du_save(:,:,id)
    407 zdvfi_omp(:,:)=dv_save(:,:,id)
    408 zdtfi_omp(:,:)=dt_save(:,:,id)
    409 pdtheta(:,:)=dtheta_save(:,:,id)
    410 zdqfi_omp(:,:,:)=dq_save(:,:,:,id)
    411 #else
    412352print*,'else'
    413353zdpsrf_omp(:)=dp_save(:)
    … …  
    415355zdvfi_omp(:,:)=dv_save(:,:)
    416356zdtfi_omp(:,:)=dt_save(:,:)
    417 pdtheta(:,:)=dtheta_save(:,:)
    418357zdqfi_omp(:,:,:)=dq_save(:,:,:)
    419 #endif
    420358!!!
    421359!!! BIG LOOP : 2. calculate physical tendencies
    422360!!!
    423 ELSE
     361ELSE ! if (test.EQ.0)
    424362!----------!
    425363! ALLOCATE !
    … …  
    489427!--------------------------------------!
    490428dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
    491 !p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> zplay_omp
     429p_prof(:) = p_hyd(i,kps:kpe,j)         ! hydrostatic pressure at layers >> zplay_omp
     430p8w_prof(:) = p_hyd_w(i,kps:kpe+1,j)         ! hydrostatic pressure at levels
     431!JL22 using hydrostatic pressures to conserve mass
    492432t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
    493433t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
    … …  
    495435v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
    496436z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> zphi_omp/g
    497 
    498 
    499 !------------------------------------------------------------!
    500 ! reconstructing hydrostatic level pressure to conserve mass !
    501 !------------------------------------------------------------!
    502 p8w_prof(kpe+1) = ptop
    503 if (TRACER_MODE .ge. 42) then
    504   do k=kpe,kps,-1
    505     p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j) * (1.d0 + moist(i,k,j,P_QV))
    506     p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
    507   enddo
    508 else
    509   do k=kpe,kps,-1
    510     p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j)
    511     p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
    512   enddo
    513 endif
    514 !p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> zplev_omp
    515437
    516438
    … …  
    644566            M_TI,CST_TI,&
    645567            M_ISOIL,M_DSOIL,&
    646             M_TSOIL,M_TSURF)
     568            M_TSOIL,P_TSURF)
    647569!!!
    648570CALL update_inputs_physiq_surf( &
    … …  
    651573            JULYR,TRACER_MODE,&
    652574            M_ALBEDO,CST_AL,&
    653             M_TSURF,M_EMISS,M_CO2ICE,&
     575            P_TSURF,M_EMISS,M_CO2ICE,&
    654576            M_GW,M_Z0,CST_Z0,&
    655577            M_H2OICE,&
    … …  
    678600!!!!!!!!!!!!!!!!!!!!!!
    679601
    680 call_physics : IF (wappel_phys .ne. 0.) THEN
    681602!!! initialize tendencies (security)
    682603zdpsrf_omp(:)=0.
    … …  
    684605zdvfi_omp(:,:)=0.
    685606zdtfi_omp(:,:)=0.
    686 pdtheta(:,:)=0.
    687607zdqfi_omp(:,:,:)=0.
     608
     609call_physics : IF (wappel_phys .ne. 0.) THEN
    688610!print *, '** ',planet_type,'** CALL TO LMD PHYSICS'
    689611!!!
    … …  
    699621!!!
    700622
    701 !! specific scenario. necessary to add the right amount of dust.
    702 #ifdef DUSTSTORM
    703 IF (firstcall .EQV. .true.) THEN
    704   zdqfi_omp(:,:,:) = zdqfi_omp(:,:,:) / dt
    705 ENDIF
    706 #endif
    707 
    708 IF (planet_type .eq. "venus" ) THEN
    709   DO j=jps,jpe
    710   DO i=ips,ipe
    711     do k=kps,kpe
    712       subs=(j-jps)*(ipe-ips+1)+(i-ips+1)
    713       tk1=(ztfi_omp(subs,k)**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
    714       tk2=((ztfi_omp(subs,k) + zdtfi_omp(subs,k))**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
    715       pdtheta(subs,k)=tk2-tk1
    716     enddo
    717   ENDDO
    718   ENDDO
    719 ENDIF
    720 
    721 !print *, '** ',planet_type,'** CALL TO LMD PHYSICS DONE'
    722 
    723623!---------------------------------------------------------------------------------!
    724624! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
    725625!---------------------------------------------------------------------------------!
    726 #ifdef SPECIAL_NEST_SAVE
    727 dp_save(:,id)=zdpsrf_omp(:)
    728 du_save(:,:,id)=zdufi_omp(:,:)
    729 dv_save(:,:,id)=zdvfi_omp(:,:)
    730 dt_save(:,:,id)=zdtfi_omp(:,:)
    731 dtheta_save(:,:,id)=pdtheta(:,:)
    732 dq_save(:,:,:,id)=zdqfi_omp(:,:,:)
    733 #else
    734626dp_save(:)=zdpsrf_omp(:)
    735627du_save(:,:)=zdufi_omp(:,:)
    736628dv_save(:,:)=zdvfi_omp(:,:)
    737629dt_save(:,:)=zdtfi_omp(:,:)
    738 dtheta_save(:,:)=pdtheta(:,:)
    739630dq_save(:,:,:)=zdqfi_omp(:,:,:)
    740 #endif
    741631
    742632!! OUTPUT OUTPUT OUTPUT
    … …  
    749639            ips,ipe,jps,jpe,&
    750640            TRACER_MODE,&
    751             M_TSURF,M_CO2ICE,&
     641            P_TSURF,M_CO2ICE,&
    752642            M_H2OICE)
    753643!!!
    … …  
    812702
    813703    ! zonal wind
    814   RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) 
     704  RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe)
    815705    ! meridional wind
    816706  RVBLTEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe)
    817707    ! potential temperature
    818708    ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
    819   IF (planet_type .eq. "venus" ) THEN
    820     RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe)
    821   ELSE
    822     RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
    823   ENDIF
     709  RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
    824710    ! update surface pressure (cf CO2 cycle in physics)
    825711    ! here dt is needed
    … …  
    841727      scalar(i,kps:kpe,j,P_QH2O_ICE)=scalar(i,kps:kpe,j,P_QH2O_ICE) &
    842728           +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
    843       ! if you want to use this mode, RTHBLTEN should be corrected as below.
    844       ! we keep it like that for the moment for testing.
     729       ! if you want to use this mode, RTHBLTEN should be corrected as below.
     730       ! we keep it like that for the moment for testing.
    845731    CASE(43)
    846732      moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
    … …  
    864750ENDDO
    865751CALL deallocate_interface
    866 DEALLOCATE(pdtheta)
    867752!!*****!!
    868753!! END !!
Note: See TracChangeset for help on using the changeset viewer.