Changeset 2868 for trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F

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.

File:

• trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F (modified) (27 diffs)

trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F

@@ -26 +26 @@
         P8W,DZ8W,T8W,Z,HT,MUT,DNW, &
         U,V,TH,T,P,EXNER,RHO, &
+        P_HYD, P_HYD_W, &
         PTOP, &
         RADT, &
@@ -35 +36 @@
         planet_type, &
         M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
-        M_H2OICE,M_TSOIL,M_Q2,M_TSURF, &
+        M_H2OICE,M_TSOIL,M_Q2,P_TSURF, &
         M_FLUXRAD,M_WSTAR,M_ISOIL,M_DSOIL,&
         M_Z0, CST_Z0, M_GW, &
@@ -41 +42 @@
         CST_AL, CST_TI, &
         isfflx, diff_opt, km_opt, &
-        HISTORY_INTERVAL, &
-        HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,&
+        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,&
@@ -94 +94 @@
 REAL, INTENT(IN  ) :: CST_AL, CST_TI 
 REAL, INTENT(IN  ) :: PTOP
-INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
 ! 2D arrays          
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
@@ -103 +102 @@
      SLPX,SLPY, &
      M_CO2ICE,M_H2OICE, &
-     M_TSURF, M_Z0, &
+     P_TSURF, M_Z0, &
      M_FLUXRAD,M_WSTAR, &
      PSFC,TSK
@@ -116 +115 @@
 ! 3D arrays 
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
-     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th
+     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,DT_AJS,RDUST,VMR_ICE,RICE,&
+     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 ) :: &
@@ -156 +155 @@
    ! <------ outputs:
    !     physical tendencies
-   REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdtheta
    ! ... intermediate arrays
    REAL, DIMENSION(:), ALLOCATABLE  :: &
@@ -171 +169 @@
    LOGICAL, SAVE :: firstcall = .true.
    INTEGER, SAVE :: previous_id = 0
-!**************************************************
-! IMPORTANT: pour travailler avec grid nesting
-! --- deux comportements distincts du save
-! ... ex1: ferme planeto avec PGF+MPI: standard
-! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE 
-!**************************************************
-#ifdef SPECIAL_NEST_SAVE
-      !  une dimension supplementaire liee au nest
-      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
-             dp_save
-      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
-             du_save, dv_save, dt_save,dtheta_save
-      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
-             dq_save      
-#else
-      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
-             dp_save
-      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
-             du_save, dv_save, dt_save,dtheta_save
-      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
-             dq_save      
-#endif
+   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: dp_save
+   REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: du_save, dv_save, dt_save
+   REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq_save      
 
 !!!IDEALIZED IDEALIZED
@@ -270 +249 @@
 ! DIMENSIONS FOR THE PHYSICS !
 !----------------------------!
-!day_ini = JULDAY - 1      !! GCM convention  !! pday at firstcall is day_ini
 wappel_phys = RADT
 zdt_split = dt*wappel_phys            ! physical timestep (s)
@@ -282 +260 @@
 lastcall = .false.
 ! **** needed but hardcoded
-
-          !PRINT *, ips, ipe, jps, jpe
-          !PRINT *, ngrid
 
 elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
@@ -300 +275 @@
   !firstcall=.false. !! just in case you'd want to get rid of the physics 
 test=0
-#ifdef SPECIAL_NEST_SAVE
-PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
-IF( .NOT. ALLOCATED( dp_save  ) ) THEN
-   PRINT *, '**** check **** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
-   !! here are the arrays that would be useful to save physics tendencies
-   ALLOCATE(dp_save(ngrid,max_dom))
-   ALLOCATE(du_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dv_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dt_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dtheta_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
-   dp_save(:,:)=0.    !! initialize these arrays ...
-   du_save(:,:,:)=0.
-   dv_save(:,:,:)=0.
-   dt_save(:,:,:)=0.
-   dtheta_save(:,:,:)=0.
-   dq_save(:,:,:,:)=0.
-ENDIF
-IF (id .lt. max_dom) THEN
-   flag_first_restart=.true.
-ELSE
-   flag_first_restart=.false.
-ENDIF
-#else
 IF( .NOT. ALLOCATED( dp_save  ) ) THEN
 ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
@@ -329 +280 @@
 ALLOCATE(dv_save(ngrid,nlayer))
 ALLOCATE(dt_save(ngrid,nlayer))
-ALLOCATE(dtheta_save(ngrid,nlayer))
 ALLOCATE(dq_save(ngrid,nlayer,nq))
 ENDIF
@@ -336 +286 @@
 dv_save(:,:)=0.
 dt_save(:,:)=0.
-dtheta_save(:,:)=0.
 dq_save(:,:,:)=0.
 flag_first_restart=.false.
-#endif
-ELSE
+
+ELSE ! is_first_step
 !-------------------------------------------------!
 ! what is done for the other steps of simulation  !
@@ -396 +345 @@
 !----------!
 CALL allocate_interface(ngrid,nlayer,nq)
-ALLOCATE(pdtheta(ngrid,nlayer))
 !!!
 !!! BIG LOOP : 1. no call for physics, used saved values
@@ -402 +350 @@
 IF (test.NE.0) THEN
 print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test
-#ifdef SPECIAL_NEST_SAVE
-zdpsrf_omp(:)=dp_save(:,id)
-zdufi_omp(:,:)=du_save(:,:,id)
-zdvfi_omp(:,:)=dv_save(:,:,id)
-zdtfi_omp(:,:)=dt_save(:,:,id)
-pdtheta(:,:)=dtheta_save(:,:,id)
-zdqfi_omp(:,:,:)=dq_save(:,:,:,id)
-#else
 print*,'else'
 zdpsrf_omp(:)=dp_save(:)
@@ -415 +355 @@
 zdvfi_omp(:,:)=dv_save(:,:)
 zdtfi_omp(:,:)=dt_save(:,:)
-pdtheta(:,:)=dtheta_save(:,:)
 zdqfi_omp(:,:,:)=dq_save(:,:,:)
-#endif
 !!!
 !!! BIG LOOP : 2. calculate physical tendencies
 !!!
-ELSE
+ELSE ! if (test.EQ.0)
 !----------!
 ! ALLOCATE !
@@ -489 +427 @@
 !--------------------------------------!
 dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
-!p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> zplay_omp
+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)
@@ -495 +435 @@
 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
-
-
-!------------------------------------------------------------!
-! reconstructing hydrostatic level pressure to conserve mass !
-!------------------------------------------------------------!
-p8w_prof(kpe+1) = ptop
-if (TRACER_MODE .ge. 42) then
-  do k=kpe,kps,-1
-    p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j) * (1.d0 + moist(i,k,j,P_QV))
-    p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
-  enddo
-else
-  do k=kpe,kps,-1
-    p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j)
-    p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
-  enddo
-endif
-!p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> zplev_omp
 
 
@@ -644 +566 @@
             M_TI,CST_TI,&
             M_ISOIL,M_DSOIL,&
-            M_TSOIL,M_TSURF)
+            M_TSOIL,P_TSURF)
 !!!
 CALL update_inputs_physiq_surf( &
@@ -651 +573 @@
             JULYR,TRACER_MODE,&
             M_ALBEDO,CST_AL,&
-            M_TSURF,M_EMISS,M_CO2ICE,&
+            P_TSURF,M_EMISS,M_CO2ICE,&
             M_GW,M_Z0,CST_Z0,&
             M_H2OICE,&
@@ -678 +600 @@
 !!!!!!!!!!!!!!!!!!!!!!
 
-call_physics : IF (wappel_phys .ne. 0.) THEN
 !!! initialize tendencies (security)
 zdpsrf_omp(:)=0.
@@ -684 +605 @@
 zdvfi_omp(:,:)=0.
 zdtfi_omp(:,:)=0.
-pdtheta(:,:)=0.
 zdqfi_omp(:,:,:)=0.
+
+call_physics : IF (wappel_phys .ne. 0.) THEN
 !print *, '** ',planet_type,'** CALL TO LMD PHYSICS'
 !!!
@@ -699 +621 @@
 !!!
 
-!! specific scenario. necessary to add the right amount of dust.
-#ifdef DUSTSTORM
-IF (firstcall .EQV. .true.) THEN
-  zdqfi_omp(:,:,:) = zdqfi_omp(:,:,:) / dt
-ENDIF
-#endif
-
-IF (planet_type .eq. "venus" ) THEN
-  DO j=jps,jpe
-  DO i=ips,ipe
-    do k=kps,kpe
-      subs=(j-jps)*(ipe-ips+1)+(i-ips+1)
-      tk1=(ztfi_omp(subs,k)**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
-      tk2=((ztfi_omp(subs,k) + zdtfi_omp(subs,k))**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
-      pdtheta(subs,k)=tk2-tk1
-    enddo
-  ENDDO
-  ENDDO
-ENDIF
-
-!print *, '** ',planet_type,'** CALL TO LMD PHYSICS DONE'
-
 !---------------------------------------------------------------------------------!
 ! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
 !---------------------------------------------------------------------------------!
-#ifdef SPECIAL_NEST_SAVE
-dp_save(:,id)=zdpsrf_omp(:)
-du_save(:,:,id)=zdufi_omp(:,:)
-dv_save(:,:,id)=zdvfi_omp(:,:)
-dt_save(:,:,id)=zdtfi_omp(:,:)
-dtheta_save(:,:,id)=pdtheta(:,:)
-dq_save(:,:,:,id)=zdqfi_omp(:,:,:)
-#else
 dp_save(:)=zdpsrf_omp(:)
 du_save(:,:)=zdufi_omp(:,:)
 dv_save(:,:)=zdvfi_omp(:,:)
 dt_save(:,:)=zdtfi_omp(:,:)
-dtheta_save(:,:)=pdtheta(:,:)
 dq_save(:,:,:)=zdqfi_omp(:,:,:)
-#endif
 
 !! OUTPUT OUTPUT OUTPUT
@@ -749 +639 @@
             ips,ipe,jps,jpe,&
             TRACER_MODE,&
-            M_TSURF,M_CO2ICE,&
+            P_TSURF,M_CO2ICE,&
             M_H2OICE)
 !!!
@@ -812 +702 @@
 
     ! zonal wind
-  RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) 
+  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)
-  IF (planet_type .eq. "venus" ) THEN
-    RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe)
-  ELSE
-    RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
-  ENDIF
+  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
@@ -841 +727 @@
       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.
+       ! 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) &
@@ -864 +750 @@
 ENDDO
 CALL deallocate_interface
-DEALLOCATE(pdtheta)
 !!*****!!
 !! END !!