source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 607

Last change on this file since 607 was 528, checked in by aslmd, 13 years ago

LMDZ.MARS : problem with r520. commits were made on an older version. fix the (numerous) problems this caused. please developers use svn update before committing.
File size: 14.4 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      IMPLICIT NONE
5c=======================================================================
6c
7c   creation des calottes pour l'etat initial
8c
9c=======================================================================
10c-----------------------------------------------------------------------
11c   Declarations:
12c   -------------
13#include "dimensions.h"
14#include "dimphys.h"
15#include "surfdat.h"
16#include "callkeys.h"
17#include "tracer.h"
18#include "comgeomfi.h"
19#include "comcstfi.h"
20c
21      INTEGER ngrid,ig,icap,iq,alternate
22      REAL  piceco2(ngrid),psolaralb(ngrid,2)
23      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
24      REAL icedryness ! ice dryness
25
26      EXTERNAL ISMIN,ISMAX
27      INTEGER ISMIN,ISMAX
28     
29! Choose false to have a somewhat non resolution dependant water ice caps distribution,
30! i.e based only on lat & lon values of each physical point.
31! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48
32! For vizualisation : soon ...
33       LOGICAL,SAVE :: improvedicelocation = .true.
34c
35c=======================================================================
36
37c     water ice outliers
38c     ------------------------------------------
39
40      IF ((water) .and. (caps)) THEN
41     
42c Perennial H20 north cap defined by watercaptag=true (allows surface to be
43c hollowed by sublimation in vdifc).
44
45c We might not want albedodat to be modified because it is used to write
46c restart files. Instead, albedo is directly modified when needed (i.e.
47c if we have watercaptag and no co2 ice), below and in albedocaps.F90
48
49c       "Dryness coefficient" controlling the evaporation and
50c        sublimation from the ground water ice (close to 1)
51c        HERE, the goal is to correct for the fact
52c        that the simulated permanent water ice polar caps
53c        is larger than the actual cap and the atmospheric
54c        opacity not always realistic.
55
56         alternate = 0
57         temptag = .true.
58         
59         write(*,*) "Ice dryness ?"
60         icedryness=1. ! default value
61         call getin("icedryness",icedryness)
62         write(*,*) " icedryness = ",icedryness
63
64       
65#ifdef MESOSCALE
66      do ig=1,ngridmx
67         !write(*,*) "all qsurf to zero. dirty."
68         do iq=1,nqmx
69         qsurf(ig,iq)=0.  !! on jette les inputs GCM
70                          !! on regle juste watercaptag
71                          !! il faudrait garder les inputs GCM
72                          !! si elles sont consequentes
73         enddo
74         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
75     .        ( albedodat(ig) .ge. 0.26   ) )  then
76                 write(*,*)"outlier ",ig
77                 watercaptag(ig)  = .true.
78                 dryness(ig)      = 1.
79                 albedodat(ig)    = albedo_h2o_ice  !! pour output
80         else
81                 watercaptag(ig)  = .false.
82                 dryness(ig)      = 1.
83         endif 
84      enddo
85#else
86       ! print*,'ngridmx',ngridmx
87
88      if (improvedicelocation) then
89     
90        if (ngridmx .eq. 738) then ! hopefully 32x24
91          print*,'water ice caps distribution for 32x24 resolution:'
92          watercaptag(1:9)    = .true. ! central cap - core
93          watercaptag(26:33)  = .true. ! central cap
94           
95        else if (ngridmx .eq. 3010) then ! hopefully 64x48
96          print*,'water ice caps distribution for 64x48 resolution:'
97          watercaptag(1:65)   = .true. ! central cap - core
98!          watercaptag(85)     = .true. ! central cap
99          watercaptag(83:85)  = .true. ! central cap - BIGGER
100          watercaptag(93:114) = .true. ! central cap
101!          watercaptag(254)    = .true. ! Korolev crater
102!          watercaptag(255)    = .true. ! Korolev crater --DECALE
103          watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG
104!--------------------------------------------------------------------         
105!          watercaptag(136)    = .true. ! outlier, lat = 78.75
106!          watercaptag(138)    = .true. ! outlier, lat = 78.75
107!          watercaptag(140)    = .true. ! outlier, lat = 78.75
108!          watercaptag(142)    = .true. ! outlier, lat = 78.75
109!          watercaptag(183)    = .true. ! outlier, lat = 78.75
110!          watercaptag(185)    = .true. ! outlier, lat = 78.75
111!          watercaptag(187)    = .true. ! outlier, lat = 78.75
112!          watercaptag(189)    = .true. ! outlier, lat = 78.75
113!          watercaptag(191)    = .true. ! outlier, lat = 78.75
114!          watercaptag(193)    = .true. ! outlier, lat = 78.75
115!--------------------------------------------------------------------         
116          watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75
117          watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75
118           
119        else if (ngridmx .ne. 1) then
120      print*,'No improved ice location for this resolution!'
121      print*,'Please set improvedicelocation to false in surfini.'
122      print*,'And check lat and lon conditions for ice caps in code.'
123      call exit(1)
124         
125        endif
126       
127        ! print caps locations
128        print*,'latitude | longitude | ig'
129        do ig=1,ngridmx
130          dryness (ig) = icedryness
131          !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs
132          !!! OU remettre comme c'était avant Aymeric et pis c'est tout
133          if (watercaptag(ig)) then
134             print*,'ice water cap', lati(ig)*180./pi,
135     .              long(ig)*180./pi, ig
136          endif
137        enddo
138
139!!------------------------ test -----------------------------
140!!-----------------------------------------------------------   
141!!      else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA ....
142!!       
143!!         do ig=1,ngridmx
144!!          dryness (ig) = 1
145!!          if (albedodat(ig) .ge. 0.35) then
146!!             watercaptag(ig)    = .true.
147!!             print*,' albedo criteria ice water cap', lati(ig)*180./pi,
148!!     .              long(ig)*180./pi, ig
149!!          endif
150!!        enddo
151!!-----------------------------------------------------------   
152!!-----------------------------------------------------------   
153     
154      else
155
156         do ig=1,ngridmx
157
158         dryness (ig) = icedryness
159         
160!!c Towards olympia planitia water caps ('relatively' low latitude ones)
161!!c---------------- proposition par AS --------------------
162!!c--------------------------------------------------------
163!!c          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
164!!c     .         ( lati(ig)*180./pi      .le.  77 ) .and.
165!!c     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
166!!c     .             ( long(ig)*180./pi .le. 120. ) )
167!!c     .           .or.
168!!c     .           ( ( long(ig)*180./pi .ge. -130. ) .and.
169!!c     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
170!!c---------------- proposition par TN --------------------
171!!c---------------- HIGHLY EXPERIMENTAL -------------------
172!!c--------------------------------------------------------     
173!!       if ( ( ( lati(ig)*180./pi .ge. 73.  ) .and. ! cap #1
174!!     .           ( lati(ig)*180./pi .le. 75.1 ) .and.
175!!     .           ( long(ig)*180./pi .ge. 95.  ) .and.
176!!     .           ( long(ig)*180./pi .le. 110. ) )
177!!     .         .or.
178!!     .         ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
179!!     .           ( lati(ig)*180./pi .le. 80.  ) .and.
180!!     .           ( long(ig)*180./pi .ge. 110. ) .and.
181!!     .           ( long(ig)*180./pi .le. 140. ) )
182!!     .         .or.
183!!     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
184!!     .           ( lati(ig)*180./pi .le. 78.  ) .and.
185!!     .           ( long(ig)*180./pi .ge. 155. ) .and.
186!!     .           ( long(ig)*180./pi .le. 180. ) )
187!!     .                .or.
188!!     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
189!!     .           ( lati(ig)*180./pi .le. 72.  ) .and.
190!!     .           ( long(ig)*180./pi .ge. 163. ) .and.
191!!     .           ( long(ig)*180./pi .le. 168. ) )
192!!     .         .or.
193!!     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
194!!     .           ( lati(ig)*180./pi .le. 78.  ) .and.
195!!     .           ( long(ig)*180./pi .ge. -160.) .and.
196!!     .           ( long(ig)*180./pi .le. -120.) ) )
197!!     .         then
198       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
199     .           ( lati(ig)*180./pi .le. 80.  ) .and.
200     .           ( long(ig)*180./pi .ge. 110. ) .and.
201     .           ( long(ig)*180./pi .le. 181. ) )
202     .         .or.
203
204     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
205     .           ( lati(ig)*180./pi .le. 76.  ) .and.
206     .           ( long(ig)*180./pi .ge. 150. ) .and.
207     .           ( long(ig)*180./pi .le. 168. ) )
208     .         .or.
209     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
210     .           ( lati(ig)*180./pi .le. 80.  ) .and.
211     .           ( long(ig)*180./pi .ge. -150.) .and.
212     .           ( long(ig)*180./pi .le. -110.) ) )
213     .         then
214     
215             if (temptag) then
216             
217               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
218                if (ngridmx.ne.1) watercaptag(ig)=.true.
219!                  print*,'ice water cap', lati(ig)*180./pi,
220!     .              long(ig)*180./pi, ig
221                  !dryness(ig) = 1.
222                  alternate = 1
223                else
224                  alternate = 0
225                endif !end if alternate = 0
226               
227             else
228             
229!               if (ngridmx.ne.1) watercaptag(ig)=.true.
230!                  write(*,*) "outliers ", lati(ig)*180./pi,
231!     .              long(ig)*180./pi
232     
233             endif ! end if temptag
234             
235           endif
236
237
238c Opposite olympia planitia water cap
239c---------------- proposition par AS --------------------
240c--------------------------------------------------------
241c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
242c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
243c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
244c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
245c---------------- proposition par TN --------------------
246c--------------------------------------------------------
247           if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
248     .            ( lati(ig)*180./pi     .le.  84 ) )
249     .          .and.
250     .          ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
251     .            ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
252!     .        ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
253!     .            ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
254!     .          ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
255!     .            ( long(ig)*180./pi .le. -70. ) ) ) ) then     !!! 64x48
256              if (ngridmx.ne.1) then
257                watercaptag(ig)=.true.
258                 print*,'ice water cap', lati(ig)*180./pi,
259     .            long(ig)*180./pi, ig
260              endif
261              !dryness(ig) = 1.
262           endif
263
264c Central cap
265c---------------- anciens reglages FF -------------------
266c--------------------------------------------------------
267
268           if (lati(ig)*180./pi.gt.84) then
269             print*,'ice water cap', lati(ig)*180./pi,
270     .         long(ig)*180./pi, ig
271             if (ngridmx.ne.1) watercaptag(ig)=.true.
272             !dryness(ig) = 1.
273c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
274c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
275c               if (ngridmx.ne.1) watercaptag(ig)=.true.
276c               dryness(ig) = 1.
277c             endif
278c             if (lati(ig)*180./pi.ge.85) then
279c               if (ngridmx.ne.1) watercaptag(ig)=.true.
280c               dryness(ig) = 1.
281c             endif
282           endif  ! (lati>80 deg)
283           
284         end do ! (ngridmx)
285         
286       endif ! of if (improvedicelocation)
287#endif     
288
289       ENDIF ! (caps & water)
290
291c ===============================================================
292c      INITIAL ALBEDO
293c ===============================================================
294
295         write(*,*)"surfini: water frost thickness",
296     s     frost_albedo_threshold
297         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
298         write(*,*)"surfini: water ice TI:", inert_h2o_ice
299
300c        To start with : Initial albedo = observed dataset
301c        -------------------------------------------------
302         DO ig=1,ngrid
303              psolaralb(ig,1)=albedodat(ig)
304              psolaralb(ig,2)=albedodat(ig)
305         END DO
306         PRINT*,'minimum albedo sans water caps',
307     s     albedodat(ISMIN(ngrid,albedodat,1))
308         PRINT*,'maximum albedo sans water caps',
309     s     albedodat(ISMAX(ngrid,albedodat,1))
310
311c        initial albedo if permanent H2O ice is present
312c        ------------------------------------------------
313         IF ((water) .and. (caps)) THEN
314           DO ig=1,ngrid
315            IF (watercaptag(ig)) THEN
316              psolaralb(ig,1) = albedo_h2o_ice
317              psolaralb(ig,2) = albedo_h2o_ice
318            ENDIF
319           END DO
320           PRINT*,'minimum albedo avec water caps',
321     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
322           PRINT*,'maximum albedo avec water caps',
323     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
324
325         ENDIF
326
327c      changing initial albedo if CO2 ice is present
328c      -------------------------------------------
329
330       DO ig=1,ngrid
331         IF (piceco2(ig) .GT. 0.) THEN
332             IF(ig.GT.ngrid/2+1) THEN
333                icap=2
334             ELSE
335                icap=1
336             ENDIF
337             psolaralb(ig,1) = albedice(icap)
338             psolaralb(ig,2) = albedice(icap)   
339         END IF
340       END DO
341
342c      changing initial albedo if water ice frost is present
343c      -------------------------------------------
344       IF (water) THEN
345          do iq=1,nqmx
346c          if there is frost and surface albedo is set to albedo_h2o_ice
347           if(noms(iq).eq."h2o_ice") then
348             do ig=1,ngrid
349              if ((piceco2(ig) .eq. 0.).and.
350     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
351                     psolaralb(ig,1) = albedo_h2o_ice
352                     psolaralb(ig,2) = albedo_h2o_ice
353c                     PRINT*,'surfini.F frost',
354c     &                  lati(ig)*180./pi, long(ig)*180./pi
355               endif
356              enddo
357           endif
358          end do
359          PRINT*,'minimum albedo avec givre et co2',
360     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
361          PRINT*,'maximum albedo avec givre et co2',
362     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
363       END IF
364         
365
366      RETURN
367      END
Note: See TracBrowser for help on using the repository browser.