Changeset 3002 for LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90

Timestamp:

Oct 3, 2017, 8:18:41 AM (7 years ago)

Author:

Ehouarn Millour

Message:

Improved slab routines.
FC

File:

• LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 (modified) (23 diffs)

LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90

@@ -1 @@
-!
+!Completed
 MODULE ocean_slab_mod
 !
@@ -40 +40 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
   !$OMP THREADPRIVATE(dt_hdiff)
+  ! heat flux convergence due to GM eddy advection
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_gm
+  !$OMP THREADPRIVATE(dt_gm)
+  ! Heat Flux correction 
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_qflux
+  !$OMP THREADPRIVATE(dt_qflux)
   ! fraction of ocean covered by sea ice (sic / (oce+sic))
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
@@ -50 +56 @@
   !$OMP THREADPRIVATE(seaice)
   ! net surface heat flux, weighted by open ocean fraction
-  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
-  !$OMP THREADPRIVATE(slab_bils)
   ! slab_bils accumulated over cpl_pas timesteps
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
@@ -108 +112 @@
    LOGICAL, PUBLIC, SAVE :: slab_hdiff
    !$OMP THREADPRIVATE(slab_hdiff)
+   LOGICAL, PUBLIC, SAVE :: slab_gm
+   !$OMP THREADPRIVATE(slab_gm)
    REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
    !$OMP THREADPRIVATE(coef_hdiff)
@@ -155 +161 @@
     ENDIF
     slabh(1)=50.
+    CALL getin_p('slab_depth',slabh(1))
     IF (nslay.GT.1) THEN 
         slabh(2)=150.
@@ -176 +183 @@
     slab_ekman=0
     CALL getin_p('slab_ekman',slab_ekman)
+! GM eddy advection (2-layers only)
+    slab_gm=.FALSE.
+    CALL getin_p('slab_gm',slab_gm)
+    IF (slab_ekman.LT.2) THEN
+       slab_gm=.FALSE.
+    ENDIF 
 ! Convective adjustment
     IF (nslay.EQ.1) THEN
@@ -205 +218 @@
          (modname,'pb allocation tslab', 1)
 
-    ALLOCATE(slab_bils(klon), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slab_bils'
-       CALL abort_physic(modname,abort_message,1)
-    ENDIF
-    slab_bils(:) = 0.0   
     ALLOCATE(bils_cum(klon), stat = error)
     IF (error /= 0) THEN
@@ -252 +259 @@
     ENDIF
 
+    ALLOCATE(dt_qflux(klon,nslay), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation dt_qflux'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+    dt_qflux(:,:) = 0.0   
+
+    IF (slab_gm) THEN !GM advection
+        ALLOCATE(dt_gm(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_gm'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_gm(:,:) = 0.0   
+    ENDIF
+
     IF (slab_ekman.GT.0) THEN ! ekman transport
         ALLOCATE(dt_ekman(klon,nslay), stat = error)
@@ -276 +299 @@
     IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
       CALL gather(zmasq,zmasq_glo)
-!$OMP MASTER  ! Only master thread
+! Master thread/process only
+!$OMP MASTER  
       IF (is_mpi_root) THEN 
           CALL ini_slab_transp(zmasq_glo)
@@ -317 +341 @@
        radsol, snow, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, qflux)
+       tsurf_new, dflux_s, dflux_l, slab_bils)
     
     USE calcul_fluxs_mod
-    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
+    USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff
     USE mod_phys_lmdz_para
 
@@ -369 +393 @@
     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
+    REAL, DIMENSION(klon), INTENT(OUT)   :: slab_bils
 
 ! Local variables
@@ -378 +402 @@
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
     ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes 
-    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
+    REAL, DIMENSION(klon) :: diff_sst, diff_siv
+    REAL, DIMENSION(klon,nslay) :: lmt_bils
     ! for surface wind stress
     REAL, DIMENSION(klon) :: u0, v0
@@ -386 +411 @@
     ! horizontal diffusion and Ekman local vars
     ! dimension = global domain (klon_glo) instead of // subdomains 
-    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
     ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
     REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
@@ -444 +469 @@
 !
 !****************************************************************************************
-    lmt_bils(:)=0.
     CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) 
     ! lmt_bils and diff_sst,siv saved by limit_slab
-    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
     ! qflux = total QFlux correction (in W/m2)
+    dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
+    IF (nslay.GT.1) THEN
+       dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
+    END IF
 
 !****************************************************************************************
@@ -469 +496 @@
         CALL gather(tslab,tslab_glo)
       ! Compute horiz transport on one process only
-!$OMP MASTER  ! Only master thread
-        IF (is_mpi_root) THEN ! Only master processus          
+        IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus          
           IF (slab_hdiff) THEN 
               dt_hdiff_glo(:,:)=0.
@@ -476 +502 @@
           IF (slab_ekman.GT.0) THEN 
               dt_ekman_glo(:,:)=0.
+          END IF
+          IF (slab_gm) THEN 
+              dt_gm_glo(:,:)=0.
           END IF
           DO i=1,cpl_pas ! time splitting for numerical stability
@@ -483 +512 @@
                   CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
                 CASE (2)
-                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
                 CASE DEFAULT
                   dt_ekman_tmp(:,:)=0.
@@ -493 +522 @@
               ENDDO
               tslab_glo=tslab_glo+dt_ekman_tmp*dtime
+              IF (slab_gm) THEN ! Gent-McWilliams eddy advection
+                dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
+                ! convert dt from K.s-1.(kg.m-2) to K.s-1   
+                DO k=1,nslay
+                  dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den) 
+                END DO
+                tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+              END IF
             ENDIF
+! GM included in Ekman_2
+!            IF (slab_gm) THEN ! Gent-McWilliams eddy advection
+!              CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp)
+!              ! convert dt_gm from K.m.s-1 to K.s-1   
+!              DO k=1,nslay
+!                dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k) 
+!              END DO
+!              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+!              dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
+!            END IF
             IF (slab_hdiff) THEN ! horizontal diffusion
               ! laplacian of slab T
@@ -509 +556 @@
             END DO
           END IF
+          IF (slab_gm) THEN
+            !dt_hdiff_glo saved in W/m2
+            dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
+          END IF
           IF (slab_ekman.GT.0) THEN 
             ! dt_ekman_glo saved in W/m2
             dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
           END IF
-        END IF ! is_mpi_root
-!$OMP END MASTER
+        END IF ! master process
 !$OMP BARRIER
         ! Send new fields back to all processes
@@ -520 +570 @@
         IF (slab_hdiff) THEN
           CALL Scatter(dt_hdiff_glo,dt_hdiff)
+        END IF
+        IF (slab_gm) THEN
+          CALL Scatter(dt_gm_glo,dt_gm)
         END IF
         IF (slab_ekman.GT.0) THEN 
@@ -533 +586 @@
 ! ***********************************
       ! Add read QFlux
-      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
+      DO k=1,nslay
+        tslab(:,k)=tslab(:,k)+dt_qflux(:,k)*cyang*dtime*cpl_pas &
+                   *slabh(1)/slabh(k) 
+      END DO
       ! Add cumulated surface fluxes
       tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
@@ -913 +969 @@
     IF (ALLOCATED(tice)) DEALLOCATE(tice)
     IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
-    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
     IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
     IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
@@ -921 +976 @@
     IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
     IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
+    IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
+    IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
 
   END SUBROUTINE ocean_slab_final