!******************************************************************************* ! 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, & RTHBLTEN,RUBLTEN,RVBLTEN, & num_3d_s,SCALAR, & num_3d_m,moist, & TRACER_MODE, & planet_type, & M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, & M_H2OICE,M_TSOIL,M_Q2,P_TSURF, & M_FLUXRAD,M_WSTAR,M_ISOIL,M_DSOIL,& M_Z0, CST_Z0, M_GW, & NUM_SOIL_LAYERS, & CST_AL, CST_TI, & isfflx, diff_opt, km_opt, & HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,& 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,SWDOWNZ,& TAU_DUST,RDUST,QSURFDUST,& MTOT,ICETOT,VMR_ICE,TAU_ICE,RICE,& 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) M_ALBEDO,M_TI,M_EMISS, & SLPX,SLPY, & M_CO2ICE,M_H2OICE, & P_TSURF, M_Z0, & M_FLUXRAD,M_WSTAR, & PSFC,TSK REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: & SWDOWNZ,& TAU_DUST,QSURFDUST,& MTOT,ICETOT,TAU_ICE,& 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 ) :: & RTHBLTEN,RUBLTEN,RVBLTEN, & HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,RDUST,VMR_ICE,RICE,& CLOUDFRAC,RH,DQICE,DQVAP,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: & M_Q2 REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT ) :: & M_TSOIL,M_ISOIL, M_DSOIL REAL, INTENT(IN ) :: CST_Z0 REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN ) :: & M_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 PRINT *,'** ',planet_type,'** TRACERS NAMES' CALL update_inputs_physiq_tracers(TRACER_MODE,nq) ! **** 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 ! !------------------------------------! print *,'** ',planet_type,' ** DOMAIN',id IF (previous_id .ne. id) THEN 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 .EQ. 0) THEN q_prof(:,1)=0.95 ELSE IF (TRACER_MODE .GE. 42) 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) = SCALAR(i,kps:kpe,j,P_QH2O_ICE) / (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 (TRACER_MODE .EQ. 20) THEN IF (firstcall .EQV. .true. .and. (.not. restart)) THEN q_prof(:,:) = 0.95 ENDIF ENDIF IF (TRACER_MODE .EQ. 32) THEN IF (firstcall .EQV. .true. .and. (.not. restart)) THEN q_prof(:,7) = 0.95 !! traceurs(7) = 'co2' ENDIF 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,& M_TI,CST_TI,& M_ISOIL,M_DSOIL,& M_TSOIL,P_TSURF) !!! CALL update_inputs_physiq_surf( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& JULYR,TRACER_MODE,& M_ALBEDO,CST_AL,& P_TSURF,M_EMISS,M_CO2ICE,& M_GW,M_Z0,CST_Z0,& M_H2OICE,& phisfi_val) !!! CALL update_inputs_physiq_turb( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,& RESTART,isles,& M_Q2,M_WSTAR) !!! CALL update_inputs_physiq_rad( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& RESTART,& M_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,M_CO2ICE,& M_H2OICE) !!! CALL update_outputs_physiq_soil( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& nsoil,& M_TSOIL) !!! CALL update_outputs_physiq_rad( & ims,ime,jms,jme,& ips,ipe,jps,jpe,& M_FLUXRAD) !!! CALL update_outputs_physiq_turb( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,kps,kpe,& M_Q2,M_WSTAR,& HFMAX,ZMAX,USTM,HFX) !!! CALL update_outputs_physiq_diag( & ims,ime,jms,jme,kms,kme,& ips,ipe,jps,jpe,kps,kpe,& SWDOWNZ,TAU_DUST,QSURFDUST,& MTOT,ICETOT,TAU_ICE,& HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,& RDUST,VMR_ICE,RICE,& 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 !! !!***************************!! RTHBLTEN(ims:ime,kms:kme,jms:jme)=0. RUBLTEN(ims:ime,kms:kme,jms:jme)=0. RVBLTEN(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 RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) ! meridional wind RVBLTEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe) ! potential temperature ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible) RTHBLTEN(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(20) !! tracer mode 20 : add a passive tracer with radioactive-like decay IF ( (i == ips) .AND. (j == jps) ) print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER.' tau_decay=60.*10. !! why not make it a namelist argument? SCALAR(i,kps:kpe,j,2) = SCALAR(i,kps:kpe,j,2)*exp(-dt/tau_decay) SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer 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)) 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+moist(i,kps:kpe,j,P_QV)) ! if you want to use this mode, RTHBLTEN should be corrected as below. ! we keep it like that for the moment for testing. 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)) 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+moist(i,kps:kpe,j,P_QV)) RTHBLTEN(i,kps:kpe,j) = RTHBLTEN(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. 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