Changeset 714 for codes/icosagcm/devel/src/dissip/dissip_gcm.f90

Timestamp:

08/03/18 16:53:37 (6 years ago)

Author:

dubos

Message:

devel : backported from trunk commits r607,r648,r649,r667,r668,r669,r706

File:

• codes/icosagcm/devel/src/dissip/dissip_gcm.f90 (modified) (9 diffs)

codes/icosagcm/devel/src/dissip/dissip_gcm.f90

@@ -35 +35 @@
 !$OMP THREADPRIVATE(rayleigh_friction_type)
 !$OMP THREADPRIVATE(rayleigh_shear)
-  REAL, SAVE    :: rayleigh_tau
+  REAL, SAVE    :: rayleigh_tau, rayleigh_limlat
 !$OMP THREADPRIVATE(rayleigh_tau)
+!$OMP THREADPRIVATE(rayleigh_limlat)
   
   REAL,SAVE    :: dtdissip
@@ -103 +104 @@
       rayleigh_friction_type=2
       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+   CASE('giant_liu_schneider') 
+      rayleigh_friction_type=99
+      IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010'
    CASE DEFAULT
       IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
@@ -117 +121 @@
          STOP
       END IF
+      IF(rayleigh_friction_type == 99) THEN
+         rayleigh_limlat=0.
+         CALL getin("rayleigh_limlat",rayleigh_limlat)
+         rayleigh_limlat = rayleigh_limlat*3.14159/180.
+      ENDIF
    END IF
 
@@ -444 +453 @@
 
     IF(mintau>0) THEN
-       itau_dissip=INT(mintau/dt)
-       dtdissip=itau_dissip*dt
+       IF (itau_dissip==0) THEN
+         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
+         itau_dissip=INT(mintau/dt)
+       ENDIF
     ELSE
-       IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
+       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
        itau_dissip=100000000
     END IF
     itau_dissip=MAX(1,itau_dissip)
-    IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
-         "itau_dissip",itau_dissip," dtdissip ",dtdissip
+    dtdissip=itau_dissip*dt
+    IF (is_master) THEN
+      PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
+      PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
+    ENDIF
 
   END SUBROUTINE init_dissip 
@@ -479 +493 @@
 
     INTEGER :: ind, shear
-    INTEGER :: l,ij
+    INTEGER :: l,ij,nn
 
 !$OMP BARRIER
@@ -515 +529 @@
 
       IF(rayleigh_friction_type>0) THEN
+       IF(rayleigh_friction_type<99) THEN
          phi=f_geopot(ind)
          ue=f_ue(ind)
@@ -524 +539 @@
             ENDDO
          END DO
+       ELSE
+         ue=f_ue(ind)
+            DO ij=ij_begin,ij_end
+              nn = ij+u_right
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                !print*, "latitude", lat_e(nn)*180./3.14159
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+              nn = ij+u_lup
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+              nn = ij+u_ldown
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+            ENDDO
+       ENDIF
       END IF
    END DO
@@ -529 +562 @@
    CALL trace_end("dissip")
 
+   CALL write_dissip_tendencies
 !$OMP BARRIER
 
     CONTAINS
+
       SUBROUTINE relax(shift_t, shift_u)
         USE dcmip_initial_conditions_test_1_2_3
@@ -559 +594 @@
        END SUBROUTINE relax
       
+       SUBROUTINE write_dissip_tendencies
+         USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
+         USE wind_mod
+         USE output_field_mod
+
+         CALL transfert_request(f_due_diss1,req_e1_vect)
+         CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
+         CALL output_field("dulon_diss1",f_buf_ulon)
+         CALL output_field("dulat_diss1",f_buf_ulat)
+!
+         CALL transfert_request(f_due_diss2,req_e1_vect)
+         CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
+         CALL output_field("dulon_diss2",f_buf_ulon)
+         CALL output_field("dulat_diss2",f_buf_ulat)
+       END SUBROUTINE write_dissip_tendencies
+
   END SUBROUTINE dissip