Changeset 5325 for branches/2015/dev_r5171_CNRS_LIM3_seaicebgc/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

Timestamp:

2015-06-01T12:00:26+02:00 (9 years ago)

Author:

aumont

Message:

various bugfixes and code changes

File:

• branches/2015/dev_r5171_CNRS_LIM3_seaicebgc/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90 (modified) (7 diffs)

branches/2015/dev_r5171_CNRS_LIM3_seaicebgc/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

@@ -28 +28 @@
 
    PUBLIC   trc_sbc   ! routine called by step.F90
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt  ! vertical profile time-step, = 2 rdttra
 
    !! * Substitutions
@@ -63 +65 @@
       !
       INTEGER  ::   ji, jj, jn                                     ! dummy loop indices
-      REAL(wp) ::   zse3t                                          ! temporary scalars
+      REAL(wp) ::   zse3t, zrtrn, zratio                           ! temporary scalars
       REAL(wp) ::   zswitch, zftra, zcd, zdtra, ztfx, ztra         ! temporary scalars
       CHARACTER (len=22) :: charout
@@ -76 +78 @@
                       CALL wrk_alloc( jpi, jpj,      zsfx   )
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrd )
+      !
+      zrtrn = 1.e-15_wp
 
       SELECT CASE( nn_ice_embd )         ! levitating or embedded sea-ice option
@@ -82 +86 @@
                                          ! (2) embedded sea-ice : salt and volume fluxes and pressure
       END SELECT
+
+      IF( ln_top_euler) THEN
+         r2dt(:) =  rdttrc(:)              ! = rdttrc (use Euler time stepping)
+      ELSE
+         IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nittrc000
+            r2dt(:) =  rdttrc(:)           ! = rdttrc (restarting with Euler time stepping)
+         ELSEIF( kt <= nittrc000 + 1 ) THEN          ! at nittrc000 or nittrc000+1
+            r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog)
+         ENDIF
+      ENDIF
+
 
       IF( kt == nittrc000 ) THEN
@@ -99 +114 @@
       ENDIF
 
-      WRITE(numout,*) ' trc_sbc initial values', trn(3,2,1,2), trb(3,2,1,2), tra(3,2,1,2)
-
       ! 0. initialization
       DO jn = 1, jptra
@@ -107 +120 @@
          !                                             ! add the trend to the general tracer trend
 
-         IF ( nn_ice_tr == -1 ) THEN  ! identical concentrations in ice and ocean (old code)
+         IF ( nn_ice_tr == -1 ) THEN  ! No tracers in sea ice (null concentration in sea ice)
 
             DO jj = 2, jpj
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zse3t = 1. / fse3t(ji,jj,1)
-                  tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zsfx(ji,jj) * r1_rau0 * trn(ji,jj,1,jn) * zse3t
+                  tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zsfx(ji,jj) * r1_rau0 * trb(ji,jj,1,jn) * zse3t
                END DO
             END DO
@@ -131 +144 @@
                   ztfx  = zftra + zswitch * zcd                ! net tracer flux (+C/D if no ice/ocean mass exchange)
    
-                  zdtra = r1_rau0 * ( ztfx + zsfx(ji,jj) * trn(ji,jj,1,jn) ) * zse3t
-                  tra(ji,jj,1,jn) = MAX( tra(ji,jj,1,jn) + zdtra, 0.) ! avoid integral ocean uptake if freezing (for iron)
+                  zdtra = r1_rau0 * ( ztfx + zsfx(ji,jj) * trb(ji,jj,1,jn) ) * zse3t
+                  IF ( zdtra < 0. ) THEN
+                     zratio = -zdtra * r2dt(1) / ( trb(ji,jj,1,jn) + zrtrn )
+                     zdtra = MIN(1.0, zratio) * zdtra ! avoid negative concentrations to arise
+                  ENDIF
+                      
+                  tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zdtra 
    
                END DO