Changeset 8438


Ignore:
Timestamp:
2017-08-15T12:21:39+02:00 (3 years ago)
Author:
dford
Message:

Initial implementation of Jozef S's code to do multivariate balancing with ERSEM.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_logchlbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90

    r7803 r8438  
    2222   USE par_trc,     ONLY: jptra          ! Tracer parameters 
    2323 
     24 
    2425   IMPLICIT NONE 
    2526   PRIVATE                    
     
    2728   PUBLIC asm_logchl_bal_ersem 
    2829 
     30   REAL(wp), PARAMETER :: tiny = 1e-20 
     31 
    2932CONTAINS 
    30  
    31    SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 
    32       &                             k_maxchlinc, logchl_bkginc, logchl_balinc ) 
    33       !!--------------------------------------------------------------------------- 
    34       !!                    ***  ROUTINE asm_logchl_bal_ersem  *** 
     33    SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 
     34      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )                      
     35 
     36                                                                         
     37                                          
     38 !! INTEGER, INTENT(in   ) :: npfts 
     39 !!REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: inc_log_chlTot 
     40 
     41      LOGICAL,  INTENT(in   )                           :: ld_logchlpftinc 
     42      INTEGER,  INTENT(in   )                           :: npfts 
     43      INTEGER,  INTENT(in   )                           :: mld_choice_bgc 
     44      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts) :: logchl_bkginc 
     45      REAL(wp), INTENT(in   )                               :: k_maxchlinc 
     46      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra)   :: logchl_balinc 
     47 
    3548      !! 
    36       !! ** Purpose :   calculate increments to ERSEM from logchl increments 
    37       !! 
    38       !! ** Method  :   convert logchl increments to chl increments 
    39       !!                split between the ERSEM PFTs 
    40       !!                spread through the mixed layer 
    41       !!                [forthcoming: calculate increments to nutrients and zooplankton] 
    42       !! 
    43       !! ** Action  :   populate logchl_balinc 
    44       !! 
    45       !! References :   forthcoming... 
    46       !!--------------------------------------------------------------------------- 
    47       !! 
    48       LOGICAL,  INTENT(in   )                               :: ld_logchlpftinc 
    49       INTEGER,  INTENT(in   )                               :: npfts 
    50       INTEGER,  INTENT(in   )                               :: mld_choice_bgc 
    51       REAL(wp), INTENT(in   )                               :: k_maxchlinc 
    52       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc 
    53       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc 
    54       !! 
    55       INTEGER                      :: ji, jj, jk 
    56       INTEGER                      :: jkmax 
    57       REAL(wp)                     :: chl_tot, chl_inc 
    58       REAL(wp), DIMENSION(jpi,jpj) :: zmld 
    59       !!--------------------------------------------------------------------------- 
    60  
    61       ! Split surface logchl incs into surface Chl1-4 incs 
    62       ! 
    63       ! In order to transform logchl incs to chl incs, need to account for the background, 
    64       ! cannot simply do 10^logchl_bkginc. Need to: 
    65       ! 1) Add logchl inc to log10(background) to get log10(analysis) 
    66       ! 2) Take 10^log10(analysis) to get analysis 
    67       ! 3) Subtract background from analysis to get chl incs 
    68       ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 
    69       ! 
    70       ! Only apply increments if all of Chl1-4 background values are > 0 
    71       ! In theory, they always will be, and if any are not that's a sign 
    72       ! that something's going wrong which the assimilation might make worse 
    73       ! 
     49  
     50 INTEGER                      :: ji, jj, jk, jp 
     51 INTEGER                      :: jkmax 
     52 REAL(wp)                     :: chl_tot, chl_inc 
     53 REAL(wp), DIMENSION(jpi,jpj) :: zmld 
     54 REAL(wp), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field  !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR 
     55 REAL(wp), DIMENSION(jptra) :: threshold 
     56  
     57 
     58 
     59         ! STEP 1: TRANSFORM ASSIMILATED LOG-FIELD INCREMENTS INTO INCREMENTS AND PRODUCE THEM FOR THE SIZE-CLASSES IF NEEDED, ALSO SAVE INCREMENTS OF LOGS IN SEPARATE VARIABLES (THIS IS IN CASE ONLY TOTAL CHLOROPHYLL WAS ASSIMILATED)  
     60 
     61  IF (lwp) WRITE(numout,*) 'Code has gotten to point started' 
     62  CALL flush(numout) 
     63 
     64 
    7465      IF ( ld_logchlpftinc ) THEN 
    7566         ! 
     
    151142      ENDIF 
    152143 
    153       ! Propagate surface Chl1-4 incs through mixed layer 
    154       ! First, choose mixed layer definition 
    155       ! 
     144       
     145       
     146    
     147 
     148 
     149 
     150  IF (lwp) WRITE(numout,*) 'Code has gotten to point log_inc transformed to inc' 
     151  CALL flush(numout) 
     152   ! STEP 2: READ CROSS_COVARIANCES 
     153 
     154  CALL read_crosscovariances(cc_chl_field, npfts) 
     155 
     156  
     157 IF (lwp) WRITE(numout,*) 'Code has gotten to point cc read' 
     158 CALL flush(numout) 
     159  !! STEP 3: CALCULATE INCREMENTS OF OTHER BIOGEOCHEMICAL VARIABLES  
     160 
     161 
     162  threshold(jp_fabm_m1+jp_fabm_n3n) = 100 
     163  threshold(jp_fabm_m1+jp_fabm_n4n) = 10 
     164  threshold(jp_fabm_m1+jp_fabm_n1p) = 15 
     165  threshold(jp_fabm_m1+jp_fabm_n5s) = 80 
     166  threshold(jp_fabm_m1+jp_fabm_z4c) = 20 
     167  threshold(jp_fabm_m1+jp_fabm_z5c) = 20 
     168  threshold(jp_fabm_m1+jp_fabm_p1c) = 25 
     169  threshold(jp_fabm_m1+jp_fabm_p2c) = 25 
     170  threshold(jp_fabm_m1+jp_fabm_p3c) = 25 
     171  threshold(jp_fabm_m1+jp_fabm_p4c) = 25 
     172  threshold(jp_fabm_m1+jp_fabm_p1n) = 10 
     173  threshold(jp_fabm_m1+jp_fabm_p2n) = 10 
     174  threshold(jp_fabm_m1+jp_fabm_p3n) = 10 
     175  threshold(jp_fabm_m1+jp_fabm_p4n) = 10 
     176  threshold(jp_fabm_m1+jp_fabm_p1p) = 10 
     177  threshold(jp_fabm_m1+jp_fabm_p2p) = 10 
     178  threshold(jp_fabm_m1+jp_fabm_p3p) = 10 
     179  threshold(jp_fabm_m1+jp_fabm_p4p) = 10 
     180 
     181  
     182 
     183  
     184  DO jj = 1, jpj 
     185         DO ji = 1, jpi 
     186                 DO jp=1,jptra 
     187 
     188                   
     189            !CALL t;ge'e;fdgfgd 
     190 
     191                         IF ( (jp /= jp_fabm_m1+jp_fabm_chl1) .AND. & 
     192                            & (jp /= jp_fabm_m1+jp_fabm_chl2) .AND. & 
     193                            & (jp /= jp_fabm_m1+jp_fabm_chl3) .AND. & 
     194                            & (jp /= jp_fabm_m1+jp_fabm_chl4)) THEN  
     195             
     196                         CALL asm_inc_biogeo_ersem(npfts, jp, trn(ji,jj,1,jp), logchl_bkginc(ji,jj,:),  cc_chl_field(ji,jj,:,:), threshold(jp), & 
     197                          & logchl_balinc(ji,jj,1,jp)) 
     198 
     199 
     200                         END IF 
     201 
     202 
     203         END DO 
     204 
     205      END DO 
     206 
     207    END DO 
     208 
     209 
     210 
     211 IF (lwp) WRITE(numout,*) 'Code has gotten to point surface logs calculated' 
     212 
     213 
     214 
     215 
     216  
     217 
     218 
     219 
     220 
     221 CALL flush(numout) 
     222       !! STEP 4: PROPAGATE EACH INCREMENT THROUGH THE MIXED LAYER 
     223 
     224 
    156225      SELECT CASE( mld_choice_bgc ) 
    157226      CASE ( 1 )                   ! Turbocline/mixing depth [W points] 
     
    179248      ! Now set MLD to bottom of a level and propagate incs equally through mixed layer 
    180249      ! 
    181       DO jj = 1, jpj 
     250 
     251       DO jj = 1, jpj 
    182252         DO ji = 1, jpi 
    183253            ! 
     
    192262            ! 
    193263            DO jk = 2, jkmax 
    194                logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) 
    195                logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) 
    196                logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) 
    197                logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) 
     264 
     265               DO jp = 1, jptra 
     266 
     267               !IF (lwp) WRITE(numout,*) logchl_balinc_ersem_chl1(ji,jj,1) 
     268               !CALL flush(numout) 
     269               logchl_balinc(ji,jj,jk,jp) = logchl_balinc(ji,jj,1,jp) 
     270              
     271               END DO 
    198272            END DO 
    199             ! 
     273 
     274   
    200275         END DO 
    201276      END DO 
    202277 
    203       ! Multivariate balancing forthcoming... 
    204  
    205    END SUBROUTINE asm_logchl_bal_ersem 
     278 
     279   IF (lwp) WRITE(numout,*) 'Code has gotten to point log transformed into inc and propagated' 
     280   CALL flush(numout) 
     281 
     282END SUBROUTINE asm_logchl_bal_ersem 
     283 
     284 
     285 
     286!CONTAINS  
     287SUBROUTINE read_crosscovariances(cc_chl_field, npfts) 
     288 
     289  !!     status = nf90_open(path = filepath, cmode = nf90_nowrite, ncid = ccfile) 
     290 
     291 INTEGER, INTENT(in  ) :: npfts 
     292 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field  !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR 
     293 INTEGER :: ncid 
     294 INTEGER :: pft 
     295 CHARACTER(LEN=256)  :: path_chl(npfts) 
     296 CHARACTER(LEN=256)  :: id_chl(npfts) 
     297 
     298 IF (npfts /= 4) THEN 
     299 
     300 IF (lwp) WRITE(numout,*) 'Wrong number of Pfts for reading crosscovs files' 
     301  CALL flush(numout) 
     302 
     303 ELSE 
     304 
     305 path_chl(1) = "crosscovs_chl1.nc" 
     306 path_chl(2) = "crosscovs_chl2.nc" 
     307 path_chl(3) = "crosscovs_chl3.nc" 
     308 path_chl(4) = "crosscovs_chl4.nc" 
     309 
     310 id_chl(1) = "cc_Chl1" 
     311 id_chl(2) = "cc_Chl2" 
     312 id_chl(3) = "cc_Chl3" 
     313 id_chl(4) = "cc_Chl4" 
     314 
     315 
     316 DO pft = 1, npfts 
     317 
     318 CALL iom_open( path_chl(pft), ncid ) 
     319 CALL iom_get( ncid, jpdom_autoglo, 'cc_NO3', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n3n)) 
     320 CALL iom_get( ncid, jpdom_autoglo, 'cc_NH4', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n4n)) 
     321 CALL iom_get( ncid, jpdom_autoglo, "cc_PO4", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n1p)) 
     322 CALL iom_get( ncid, jpdom_autoglo, "cc_SiO", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n5s)) 
     323 CALL iom_get( ncid, jpdom_autoglo, "cc_Z4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z4c)) 
     324 CALL iom_get( ncid, jpdom_autoglo, "cc_Z5c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z5c)) 
     325 CALL iom_get( ncid, jpdom_autoglo, "cc_P1c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1c)) 
     326 CALL iom_get( ncid, jpdom_autoglo, "cc_P2c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2c)) 
     327 CALL iom_get( ncid, jpdom_autoglo, "cc_P3c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3c)) 
     328 CALL iom_get( ncid, jpdom_autoglo, "cc_P4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4c)) 
     329 CALL iom_get( ncid, jpdom_autoglo, "cc_P1n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1n)) 
     330 CALL iom_get( ncid, jpdom_autoglo, "cc_P2n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2n)) 
     331 CALL iom_get( ncid, jpdom_autoglo, "cc_P3n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3n)) 
     332 CALL iom_get( ncid, jpdom_autoglo, "cc_P4n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4n)) 
     333 CALL iom_get( ncid, jpdom_autoglo, "cc_P1p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1p)) 
     334 CALL iom_get( ncid, jpdom_autoglo, "cc_P2p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2p)) 
     335 CALL iom_get( ncid, jpdom_autoglo, "cc_P3p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3p)) 
     336 CALL iom_get( ncid, jpdom_autoglo, "cc_P4p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4p)) 
     337 CALL iom_get( ncid, jpdom_autoglo, id_chl(pft), cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_chl1)) 
     338 CALL iom_close( ncid ) 
     339 
     340 END DO 
     341 
     342 END IF 
     343 
     344 
     345  END SUBROUTINE read_crosscovariances 
     346 
     347 
     348!CONTAINS 
     349 SUBROUTINE asm_inc_biogeo_ersem(npfts, jp, bgf_value, chl, cc, threshold, output) 
     350 INTEGER, INTENT(in ) :: npfts 
     351 INTEGER, INTENT(in ) :: jp 
     352 REAL, INTENT(in  )  :: bgf_value 
     353 REAL, INTENT(in   ), DIMENSION(npfts) :: chl 
     354 REAL, INTENT(in   ), DIMENSION(npfts,jptra) :: cc  !! VARIANCE OF CHL SIZE-CLASS  
     355 REAL, INTENT(in   ) :: threshold 
     356 REAL, INTENT(  out) :: output 
     357 
     358 
     359 REAL :: total, weight_1, weight_2, weight_3, weight_4 
     360 
     361 
     362    output = 0.0 
     363 
     364    IF ( cc(1,jp)*cc(2,jp)*cc(3,jp)*cc(4,jp)*cc(1, jp_fabm_m1+jp_fabm_chl1)*cc(2, jp_fabm_m1+jp_fabm_chl2)*cc(3, jp_fabm_m1+jp_fabm_chl3)*cc(4, jp_fabm_m1+jp_fabm_chl4) /= 0.0  ) THEN 
     365 
     366       total = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1)) + ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2)) + ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3)) + ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4)) 
     367       weight_1 = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))/total 
     368       weight_2 = ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))/total 
     369       weight_3 = ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))/total 
     370       weight_4 = ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))/total 
     371 
     372!       PRINT *, 'weights', weight_1, weight_2, weight_3, weight_4  
     373       
     374       output = weight_1*(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))*chl(1) + weight_2*(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))*chl(2) + weight_3*(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))*chl(3) + weight_4*(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))*chl(4) 
     375 
     376 
     377       ! FIRST BREAK = DON'T INCREASE VALUE MORE THAN TENFOLD  
     378        
     379       output = MIN(output, LOG10(11.0)) 
     380       
     381 
     382       ! SECOND BREAK = IF VALUE LARGER THAN THRESHOLD DO NOT INCREASE IT 
     383 
     384       IF (bgf_value > threshold) THEN 
     385           output = MIN(output, 0.0) 
     386       ENDIF 
     387 
     388        ! THIRD BREAK = IF VALUE MORE THAN FIFTH OF THRESHOLD THAN DON'T INCREASE IT BY MORE THAN AN INCREMENT EQUAL TO TENTH OF THRESHOLD (INCREASE IT SLOWLY) 
     389 
     390       IF (bgf_value > 0.2*threshold) THEN 
     391           output = MIN(output, LOG10(0.1*threshold/bgf_value + 1)) 
     392       ENDIF 
     393 
     394       output = (10**( output ) - 1)*bgf_value 
     395   
     396       PRINT *, 'output', output 
     397          
     398       
     399    ENDIF 
     400 
     401       
     402     !! calculate increments for one or all variables 
     403 
     404   END SUBROUTINE asm_inc_biogeo_ersem 
     405 
     406 
     407 
     408   
     409 
    206410 
    207411#else 
     
    210414   !!---------------------------------------------------------------------- 
    211415CONTAINS 
     416   CONTAINS 
    212417   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 
    213418      &                             k_maxchlinc, logchl_bkginc, logchl_balinc ) 
     
    220425      WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?' 
    221426   END SUBROUTINE asm_logchl_bal_ersem 
     427 
     428 
     429 
     430 
     431 
     432   
     433 
    222434#endif 
    223  
    224435   !!====================================================================== 
    225436END MODULE asmlogchlbal_ersem 
Note: See TracChangeset for help on using the changeset viewer.