!******************************************************************************* ! PURPOSE: LMD_driver is the WRF mediation layer routine that provides the ! interface to LMD physics packages in the WRF model layer. For those familiar ! with the LMD GCM, the aim of this driver is to do part of the job of the ! calfis.F routine. !******************************************************************************* ! AUTHOR: A. Spiga - January 2007 !******************************************************************************* ! UPDATES: - included all final updates for the paper - March 2008 ! - general cleaning of code and comments - October 2008 ! - additions for idealized cases - January 2009 ! - additions for new soil model in physics - January 2010 ! - unified module_lmd_driver: old, new phys and LES - February 2011 ! - a new paradigm to prepare nested simulations - April 2014 ! - adapted to new interface philosophy & other planets - August 2016 ! - adapated to WRFV4 - JL22 !******************************************************************************* MODULE module_lmd_driver CONTAINS SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, & IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, & i_start,i_end,j_start,j_end,kts,kte,num_tiles, & DX,DY, & MSFT,MSFU,MSFV, & GMT,JULYR,JULDAY, & P8W,DZ8W,T8W,Z,HT,MUT,DNW, & U,V,TH,T,P,EXNER,RHO, & P_HYD, P_HYD_W, & PTOP, & RADT, & TSK,PSFC, & RTHPLATEN,RUPLATEN,RVPLATEN, & num_3d_s,SCALAR, & num_3d_m,moist, & TRACER_MODE, & planet_type, & P_ALBEDO,P_TI,P_CO2ICE,P_EMISS, & P_H2OICE,P_TSOIL,P_Q2,P_TSURF, & P_FLUXRAD,P_WSTAR,P_ISOIL,P_DSOIL,& P_Z0, CST_Z0, P_GW, & NUM_SOIL_LAYERS, & CST_AL, CST_TI, & isfflx, diff_opt, km_opt, & HR_SW,HR_LW,HR_DYN,DT_RAD,& CLOUDFRAC,TOTCLOUDFRAC,RH, & DQICE,DQVAP,REEVAP,SURFRAIN,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,& FLUXSURF_LW,FLXGRD,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF,LATENT_HF,& HFMAX,ZMAX,& USTM,HFX,& SLPX,SLPY,RESTART) ! NB: module_lmd_driver_output1.inc : output arguments generated from Registry !================================================================== ! USES !================================================================== USE module_state_description USE module_model_constants USE module_wrf_error !!!!!!!! interface modules USE variables_mod !! to set variables USE update_inputs_physiq_mod !! to set inputs for physiq USE update_outputs_physiq_mod !! to get outputs from physiq USE iniphysiq_mod !! to get iniphysiq subroutine USE callphysiq_mod !! to call the LMD physics !!!!!!!! interface modules !================================================================== IMPLICIT NONE !================================================================== !================================================================== ! VARIABLES !================================================================== CHARACTER(len=10),INTENT(IN) :: planet_type ! WRF Dimensions INTEGER, INTENT(IN ) :: & ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & kts,kte,num_tiles, & NUM_SOIL_LAYERS INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & i_start,i_end,j_start,j_end ! Scalars INTEGER, INTENT(IN ) :: JULDAY, itimestep, julyr,id,max_dom,TRACER_MODE INTEGER, INTENT(IN ) :: isfflx,diff_opt,km_opt REAL, INTENT(IN ) :: GMT,dt,dx,dy,RADT REAL, INTENT(IN ) :: CST_AL, CST_TI REAL, INTENT(IN ) :: PTOP ! 2D arrays REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: & MSFT,MSFU,MSFV, & XLAT,XLONG,HT, & MUT, & !total dry air column mass (in Pa) P_ALBEDO,P_TI,P_EMISS, & SLPX,SLPY, & P_CO2ICE,P_H2OICE, & P_TSURF, P_Z0, & P_FLUXRAD,P_WSTAR, & PSFC,TSK REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: & HFMAX,ZMAX,& USTM,HFX,TOTCLOUDFRAC,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,& FLUXSURF_LW,FLXGRD,LATENT_HF,REEVAP,SURFRAIN REAL, DIMENSION(kms:kme), INTENT(IN ) :: DNW ! del(eta), depth between full levels in eta variables. ! 3D arrays REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: & dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th,p_hyd,p_hyd_w REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: & RTHPLATEN,RUPLATEN,RVPLATEN, & HR_SW,HR_LW,HR_DYN,DT_RAD,& CLOUDFRAC,RH,DQICE,DQVAP,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: & P_Q2 REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT ) :: & P_TSOIL,P_ISOIL, P_DSOIL REAL, INTENT(IN ) :: CST_Z0 REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN ) :: & P_GW ! 4D arrays INTEGER, INTENT(IN ) :: num_3d_s,num_3d_m REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: scalar REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_m ), INTENT(INOUT ) :: moist ! Logical LOGICAL, INTENT(IN ) :: restart !------------------------------------------- ! LOCAL VARIABLES !------------------------------------------- INTEGER :: i,j,k,its,ite,jts,jte,ij,kk INTEGER :: subs,iii ! *** for tracer Mode 20 *** REAL :: tau_decay ! *** ----------------------- *** ! *** for LMD physics ! ------> inputs: INTEGER :: ngrid,nlayer,nq,nsoil REAL :: MY REAL :: phisfi_val LOGICAL :: lastcall ! ---------- ! <------ outputs: ! physical tendencies ! ... intermediate arrays REAL, DIMENSION(:), ALLOCATABLE :: & dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, & u_prof,v_prof,z_prof REAL, DIMENSION(:,:), ALLOCATABLE :: q_prof ! Additional control variables INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe REAL :: elaps INTEGER :: test REAL :: wappel_phys LOGICAL, SAVE :: flag_first_restart LOGICAL, SAVE :: firstcall = .true. INTEGER, SAVE :: previous_id = 0 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: dp_save REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: du_save, dv_save, dt_save REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq_save !!!IDEALIZED IDEALIZED REAL :: lat_input, lon_input, ls_input, lct_input LOGICAL :: isles !!!IDEALIZED IDEALIZED REAL :: tk1,tk2 !================================================================== ! CODE !================================================================== print *, '** ',planet_type,'** ENTER WRF-LMD DRIVER' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! determine here if this is turbulence-resolving or not IF (JULYR .le. 8999) THEN isles = .false. ! "True" LES is not available in this version IF (firstcall .EQV. .true.) PRINT *, '*** REAL-CASE SIMULATION ***' ELSE IF (firstcall .EQV. .true.) PRINT *, '*** IDEALIZED SIMULATION ***' IF ((diff_opt .eq. 2) .and. (km_opt .eq. 2)) THEN IF (firstcall .EQV. .true.) THEN PRINT *, '*** type: LES ***' PRINT *, '*** diff_opt = 2 *** km_opt = 2' PRINT *, '*** forcing is isfflx = ',isfflx ENDIF isles = .true. !! SPECIAL LES ELSE IF (firstcall .EQV. .true.) THEN PRINT *, '*** type: not LES ***' PRINT *, '*** diff_opt = ',diff_opt PRINT *, '*** km_opt = ',km_opt ENDIF isles = .false. !! IDEALIZED, no LES !! cependant, ne veut-on pas pouvoir !! prescrire un flux en idealise ?? ENDIF ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------! ! TWEAK on WRF DIMENSIONS ! !-------------------------! its = i_start(1) ! define here one big tile ite = i_end(num_tiles) jts = j_start(1) jte = j_end(num_tiles) !! IF (isles .eqv. .false.) THEN relax=0 sponge_top=0 ! another value than 0 triggers instabilities IF (id .gt. 1) relax=2 ! fix to avoid noise in nesting simulations ; 1 >> too much noise ... ENDIF ips=its ipe=ite jps=jts jpe=jte IF (isles .eqv. .false.) THEN IF (ips .eq. ids) ips=its+relax !! IF tests necesary for parallel runs IF (ipe .eq. ide-1) ipe=ite-relax IF (jps .eq. jds) jps=jts+relax IF (jpe .eq. jde-1) jpe=jte-relax ENDIF kps=kts !! start at surface IF (isles .eqv. .false.) THEN kpe=kte-sponge_top ELSE IF (firstcall .EQV. .true.) PRINT *, '*** IDEALIZED SIMULATION: LES *** kpe=kte' kpe=kte !-sponge_top ENDIF !----------------------------! ! DIMENSIONS FOR THE PHYSICS ! !----------------------------! wappel_phys = RADT zdt_split = dt*wappel_phys ! physical timestep (s) ngrid=(ipe-ips+1)*(jpe-jps+1) ! size of the horizontal grid nlayer = kpe-kps+1 ! number of vertical layers: nlayermx nsoil = NUM_SOIL_LAYERS ! number of soil layers: nsoilmx CALL update_inputs_physiq_tracers(TRACER_MODE,nq) IF (firstcall .EQV. .true.) PRINT *,'** ',planet_type,'** TRACER MODE', TRACER_MODE ! **** needed but hardcoded lastcall = .false. ! **** needed but hardcoded elaps = (float(itimestep)-1.)*dt ! elapsed seconds of simulation !----------------------------------------------! ! what is done at the first step of simulation ! !----------------------------------------------! IF (elaps .eq. 0.) THEN flag_first_restart = .false. ELSE flag_first_restart=flag_first_restart.OR.(.NOT. ALLOCATED(dp_save)) ENDIF is_first_step: IF ((elaps .eq. 0.).or.flag_first_restart) THEN firstcall=.true. !! for continuity with GCM, physics are always called at the first WRF timestep !firstcall=.false. !! just in case you'd want to get rid of the physics test=0 IF( .NOT. ALLOCATED( dp_save ) ) THEN ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies ALLOCATE(du_save(ngrid,nlayer)) ALLOCATE(dv_save(ngrid,nlayer)) ALLOCATE(dt_save(ngrid,nlayer)) ALLOCATE(dq_save(ngrid,nlayer,nq)) ENDIF dp_save(:)=0. !! initialize these arrays ... du_save(:,:)=0. dv_save(:,:)=0. dt_save(:,:)=0. dq_save(:,:,:)=0. flag_first_restart=.false. ELSE ! is_first_step !-------------------------------------------------! ! what is done for the other steps of simulation ! !-------------------------------------------------! IF (wappel_phys .gt. 0.) THEN firstcall = .false. test = MODULO(itimestep-1,int(wappel_phys)) ! WRF time is integrated time (itimestep=1 at the end of first step) ! -- same strategy as diagfi in the LMD GCM ! -- arguments of modulo must be of the same kind (here: integers) ELSE firstcall = .false. test = 9999 ENDIF ENDIF is_first_step !------------------------------------! ! print info about domain initially ! ! ... and whenever domain is changed ! !------------------------------------! IF (previous_id .ne. id) THEN print *,'** ',planet_type,' ** DOMAIN',id print *, '** ',planet_type,' ** ... INITIALIZE DOMAIN',id print *, '** ',planet_type,' ** ... PREVIOUS DOMAIN was',previous_id print *, 'ITIMESTEP: ',itimestep print *, 'TILES: ', i_start,i_end, j_start, j_end ! numbers for simple runs, arrays for parallel runs print *, 'DOMAIN: ', ids, ide, jds, jde, kds, kde print *, 'MEMORY: ', ims, ime, jms, jme, kms, kme print *, 'COMPUT: ', ips, ipe, jps, jpe, kps, kpe print *, 'SIZE OF PHYSICS GRID for this process: ',ngrid print *, 'ADVECTED TRACERS: ', num_3d_s-1 print *, 'PHYSICS IS CALLED EACH...',wappel_phys print *, '-- means: PHYSICAL TIMESTEP=',zdt_split ENDIF !!!***********!! !!! IDEALIZED !! !!!***********!! IF (.not.(JULYR .le. 8999)) THEN IF (firstcall .EQV. .true.) THEN PRINT *,'** ',planet_type,'** IDEALIZED SIMULATION' PRINT *,'** ',planet_type,'** BEWARE: input_coord must be here' ENDIF open(unit=14,file='input_coord',form='formatted',status='old') rewind(14) read(14,*) lon_input read(14,*) lat_input read(14,*) ls_input read(14,*) lct_input close(14) ENDIF !----------! ! ALLOCATE ! !----------! CALL allocate_interface(ngrid,nlayer,nq) !!! !!! BIG LOOP : 1. no call for physics, used saved values !!! IF (test.NE.0) THEN print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test print*,'else' zdpsrf_omp(:)=dp_save(:) zdufi_omp(:,:)=du_save(:,:) zdvfi_omp(:,:)=dv_save(:,:) zdtfi_omp(:,:)=dt_save(:,:) zdqfi_omp(:,:,:)=dq_save(:,:,:) !!! !!! BIG LOOP : 2. calculate physical tendencies !!! ELSE ! if (test.EQ.0) !----------! ! ALLOCATE ! !----------! ! interm ALLOCATE(dz8w_prof(nlayer)) ALLOCATE(p8w_prof(nlayer+1)) ALLOCATE(p_prof(nlayer)) ALLOCATE(t_prof(nlayer)) ALLOCATE(t8w_prof(nlayer)) ALLOCATE(u_prof(nlayer)) ALLOCATE(v_prof(nlayer)) ALLOCATE(z_prof(nlayer)) ALLOCATE(q_prof(nlayer,nq)) !!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !! PREPARE PHYSICS INPUTS !! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !! INITIALIZE AND ALLOCATE EVERYTHING !! here, only firstcall allocation_firstcall: IF (firstcall .EQV. .true.) THEN !! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM) !! parameters are defined in the module_model_constants.F WRF routine PRINT *,'** ',planet_type,'** INITIALIZE ARRAYS FOR PHYSICS' !! need to get initial time first CALL update_inputs_physiq_time(& JULYR,JULDAY,GMT,& elaps,& lct_input,lon_input,ls_input,& MY) !! Set up initial time phour_ini = JH_cur_split !! Fill planetary parameters in modules !! Values defined in the module_model_constants.F WRF routine CALL update_inputs_physiq_constants !! JL21 it seems nothing is done in update_inputs_physiq_constants for generic case !!! !! Initialize physics CALL iniphysiq(ngrid,nlayer,nq,wappel_phys,& wdaysec,floor(JD_cur), & 1./reradius,g,r_d,cp,1) ENDIF allocation_firstcall !!*****************************!! !! PREPARE "FOOD" FOR PHYSIQ.F !! !!*****************************!! DO j = jps,jpe DO i = ips,ipe !!*******************************************!! !! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !! !!*******************************************!! !-----------------------------------! ! 1D subscript for physics "cursor" ! !-----------------------------------! subs = (j-jps)*(ipe-ips+1)+(i-ips+1) !--------------------------------------! ! 1D columns : inputs for the physics ! ! ... vert. coord., meteor. fields ... ! !--------------------------------------! dz8w_prof(:) = dz8w(i,kps:kpe,j) ! dz between full levels (m) p_prof(:) = p_hyd(i,kps:kpe,j) ! hydrostatic pressure at layers >> zplay_omp p8w_prof(:) = p_hyd_w(i,kps:kpe+1,j) ! hydrostatic pressure at levels !JL22 using hydrostatic pressures to conserve mass t_prof(:) = t(i,kps:kpe,j) ! temperature half level (K) >> pt t8w_prof(:) = t8w(i,kps:kpe,j) ! temperature full level (K) u_prof(:) = u(i,kps:kpe,j) ! zonal wind (A-grid: unstaggered) half level (m/s) >> zufi_omp v_prof(:) = v(i,kps:kpe,j) ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv z_prof(:) = z(i,kps:kpe,j) ! geopotential height half level (m) >> zphi_omp/g !--------------------------------! ! specific treatment for tracers ! !--------------------------------! IF (TRACER_MODE == 1) THEN ! to be clean we should have an automatized process that makes sure that moist is sent to igcm_h2o_vap and etc. q_prof(:,1) = SCALAR(i,kps:kpe,j,P_QH2O) / (1.d0 + SCALAR(i,kps:kpe,j,P_QH2O)) !! P_xxx is the index for variable xxx. q_prof(:,2) = SCALAR(i,kps:kpe,j,P_QH2O_ICE) / (1.d0 + SCALAR(i,kps:kpe,j,P_QH2O)) ! conversion from mass mixing ratio in WRF to specific concentration in Physiq ELSE IF ((TRACER_MODE >= 42).AND.(TRACER_MODE <= 45)) THEN ! to be clean we should have an automatized process that makes sure that moist is sent to igcm_h2o_vap and etc. q_prof(:,1) = moist(i,kps:kpe,j,P_QV) / (1.d0 + moist(i,kps:kpe,j,P_QV)) !! P_xxx is the index for variable xxx. q_prof(:,2) = moist(i,kps:kpe,j,P_QC) / (1.d0 + moist(i,kps:kpe,j,P_QV)) ! conversion from mass mixing ratio in WRF to specific concentration in Physiq ELSE q_prof(:,1:nq) = SCALAR(i,kps:kpe,j,2:nq+1) !! the names were set above !! one dummy tracer in WRF !JL22 cannot normalize to moist here as we do not know if it has been initialized. ENDIF IF (firstcall .EQV. .true.) THEN !-----------------! ! Optional output ! !-----------------! IF ( (i == ips) .AND. (j == jps) ) THEN PRINT *,'z_prof ',z_prof PRINT *,'dz8w_prof',dz8w_prof PRINT *,'p8w_prof ',p8w_prof PRINT *,'p_prof ',p_prof PRINT *,'t_prof ',t_prof PRINT *,'t8w_prof ',t8w_prof PRINT *,'u_prof ',u_prof PRINT *,'v_prof ',v_prof PRINT *,'q_prof ',q_prof ENDIF ENDIF !---------------------------------------------! ! in LMD physics, geopotential must be ! ! expressed with respect to the local surface ! !---------------------------------------------! zphi_omp(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) ) !--------------------------------! ! Dynamic fields for LMD physics ! !--------------------------------! zplev_omp(subs,1:nlayer+1) = p8w_prof(1:nlayer+1) !! NB: last level: no data zplay_omp(subs,:) = p_prof(:) ztfi_omp(subs,:) = t_prof(:) zufi_omp(subs,:) = u_prof(:) zvfi_omp(subs,:) = v_prof(:) flxwfi_omp(subs,:) = 0 !! NB: not used in the physics, only diagnostic... !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! for IDEALIZED CASES ONLY !IF (.not.(JULYR .le. 8999)) zplev_omp(subs,nlayer+1)=0. !! zplev_omp(subs,nlayer+1)=ptop >> NO ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! NOTE: ! IF ( zplev_omp(subs,nlayer+1) .le. 0 ) zplev_omp(subs,nlayer+1)=ptop ! cree des diagnostics delirants et aleatoires dans le transfert radiatif !---------! ! Tracers ! !---------! zqfi_omp(subs,:,:) = q_prof(:,:) !! traceurs generiques, seuls noms sont specifiques ENDDO ENDDO !---------------------------------------------------------! ! Ground geopotential ! ! (=g*HT(i,j), only used in the microphysics: newcondens) ! !---------------------------------------------------------! phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.) !! NB: z_prof are half levels, so z_prof(1) is not the surface !!*****************!! !! CLEAN THE PLACE !! !!*****************!! DEALLOCATE(dz8w_prof) DEALLOCATE(z_prof) DEALLOCATE(p8w_prof) DEALLOCATE(p_prof) DEALLOCATE(t_prof) DEALLOCATE(t8w_prof) DEALLOCATE(u_prof) DEALLOCATE(v_prof) DEALLOCATE(q_prof) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! ONE DOMAIN CASE !!! --> we simply need to initialize static and physics inputs !!! SEVERAL DOMAINS !!! --> we update static and physics inputs each time !!! to account for domain change !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN PRINT *,'** ',planet_type,'** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.' !!*******************************************!! !!*******************************************!! CALL update_inputs_physiq_geom( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& JULYR,ngrid,nlayer,& DX,DY,MSFT,& lat_input, lon_input,& XLAT,XLONG) !!! CALL update_inputs_physiq_slope( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& JULYR,& SLPX,SLPY) !!! print*, 'num_soil_layers, nsoil', num_soil_layers, nsoil CALL update_inputs_physiq_soil( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& JULYR,nsoil,& P_TI,CST_TI,& P_ISOIL,P_DSOIL,& P_TSOIL,P_TSURF) !!! CALL update_inputs_physiq_surf( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& JULYR,TRACER_MODE,& P_ALBEDO,CST_AL,& P_TSURF,P_EMISS,P_CO2ICE,& P_GW,P_Z0,CST_Z0,& P_H2OICE,& phisfi_val) !!! CALL update_inputs_physiq_turb( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,& RESTART,isles,& P_Q2,P_WSTAR) !!! CALL update_inputs_physiq_rad( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& RESTART,& P_FLUXRAD) !!! ENDIF pass_interface !!*******************************************!! !!*******************************************!! !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! !! CALL LMD PHYSICS !! !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! !!! initialize tendencies (security) zdpsrf_omp(:)=0. zdufi_omp(:,:)=0. zdvfi_omp(:,:)=0. zdtfi_omp(:,:)=0. zdqfi_omp(:,:,:)=0. call_physics : IF (wappel_phys .ne. 0.) THEN !print *, '** ',planet_type,'** CALL TO LMD PHYSICS' !!! CALL update_inputs_physiq_time(& JULYR,JULDAY,GMT,& elaps,& lct_input,lon_input,ls_input,& MY) !!! CALL call_physiq(planet_type,ngrid,nlayer,nq, & firstcall,lastcall) !!! !---------------------------------------------------------------------------------! ! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS ! !---------------------------------------------------------------------------------! dp_save(:)=zdpsrf_omp(:) du_save(:,:)=zdufi_omp(:,:) dv_save(:,:)=zdvfi_omp(:,:) dt_save(:,:)=zdtfi_omp(:,:) dq_save(:,:,:)=zdqfi_omp(:,:,:) !! OUTPUT OUTPUT OUTPUT !-------------------------------------------------------! ! Save key variables for restart and output and nesting ! !-------------------------------------------------------! !!! CALL update_outputs_physiq_surf( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& TRACER_MODE,& P_TSURF,P_CO2ICE,& P_H2OICE) !!! CALL update_outputs_physiq_soil( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& nsoil,& P_TSOIL) !!! CALL update_outputs_physiq_rad( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& P_FLUXRAD) !!! CALL update_outputs_physiq_turb( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,kps,kpe,& P_Q2,P_WSTAR,& HFMAX,ZMAX,USTM,HFX) !!! CALL update_outputs_physiq_diag( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,kps,kpe,& HR_SW,HR_LW,HR_DYN,DT_RAD,& CLOUDFRAC,TOTCLOUDFRAC,RH,& DQICE,DQVAP,REEVAP,SURFRAIN,& ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,& FLUXSURF_LW,FLXGRD,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF,LATENT_HF) !!! !print *, '** ',planet_type,'** OUTPUT PHYSICS DONE' ENDIF call_physics ENDIF ! end of BIG LOOP 2. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!***************************!! !! DEDUCE TENDENCIES FOR WRF !! !!***************************!! RTHPLATEN(ims:ime,kms:kme,jms:jme)=0. RUPLATEN(ims:ime,kms:kme,jms:jme)=0. RVPLATEN(ims:ime,kms:kme,jms:jme)=0. PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF !------------------------------------------------------------------! ! WRF adds the physical and dynamical tendencies ! ! thus the tendencies must be passed as 'per second' tendencies ! ! --when physics is not called, the applied physical tendency ! ! --is the one calculated during the last call to physics ! !------------------------------------------------------------------! !print*,'pdt',pdt(1,1),pdt(1,nlayer) !print*,'exner',exner(1,:,1) DO j = jps,jpe DO i = ips,ipe subs = (j-jps)*(ipe-ips+1)+(i-ips+1) ! zonal wind RUPLATEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) ! meridional wind RVPLATEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe) ! potential temperature ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible) RTHPLATEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j) ! update surface pressure (cf CO2 cycle in physics) ! here dt is needed PSFC(i,j)=PSFC(i,j)+zdpsrf_omp(subs)*dt ! tracers SCALAR(i,kps:kpe,j,1)=0. SELECT CASE (TRACER_MODE) CASE(0) SCALAR(i,kps:kpe,j,:)=0. CASE(1) scalar(i,kps:kpe,j,P_QH2O)=scalar(i,kps:kpe,j,P_QH2O) & +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+scalar(i,kps:kpe,j,P_QH2O)) scalar(i,kps:kpe,j,P_QH2O_ICE)=scalar(i,kps:kpe,j,P_QH2O_ICE) & +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+scalar(i,kps:kpe,j,P_QH2O)) ! if you want to use this mode, RTHPLATEN should be corrected as below. ! we keep it like that for the moment for testing. CASE(42) moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) & +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) & +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) RTHPLATEN(i,kps:kpe,j) = RTHPLATEN(i,kps:kpe,j) & * (1.d0+moist(i,kps:kpe,j,P_QV))/(1.d0+rvovrd*moist(i,kps:kpe,j,P_QV)) ! correct dT/dt assuming a constant molar heat capacity. ! Specific heat cappacity scales with molar mass. CASE(43) moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) & +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) & +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) tau_decay=86400.*100. !! why not make it a namelist argument? SCALAR(i,kps:kpe,j,P_MARKER) = SCALAR(i,kps:kpe,j,P_MARKER)*exp(-dt/tau_decay) SCALAR(i,1,j,P_MARKER) = 1. !! this tracer is emitted in the surface layer RTHPLATEN(i,kps:kpe,j) = RTHPLATEN(i,kps:kpe,j) & * (1.d0+moist(i,kps:kpe,j,P_QV))/(1.d0+rvovrd*moist(i,kps:kpe,j,P_QV)) ! correct dT/dt assuming a constant molar heat capacity. ! Specific heat cappacity scales with molar mass. CASE(44) moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) & +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) & +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) CASE(45) moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) & +zdqfi_omp(subs,kps:kpe,1)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) moist(i,kps:kpe,j,P_QC)=moist(i,kps:kpe,j,P_QC) & +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV)) tau_decay=86400.*100. !! why not make it a namelist argument? SCALAR(i,kps:kpe,j,P_MARKER) = SCALAR(i,kps:kpe,j,P_MARKER)*exp(-dt/tau_decay) SCALAR(i,1,j,P_MARKER) = 1. !! this tracer is emitted in the surface layer CASE DEFAULT !SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+zdqfi_omp(subs,kps:kpe,1:nq)*dt !!! here dt is needed scalar(i,kps:kpe,j,2:nq+1)=scalar(i,kps:kpe,j,2:nq+1) & +zdqfi_omp(subs,kps:kpe,1:nq)*dt END SELECT ENDDO ENDDO CALL deallocate_interface !!*****!! !! END !! !!*****!! !print *,'** ',planet_type,'** END LMD PHYSICS' previous_id = id !----------------------------------------------------------------! ! use debug (see solve_em) if you wanna display some message ... ! !----------------------------------------------------------------! END SUBROUTINE lmd_driver END MODULE module_lmd_driver