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

Last change on this file since 520 was 520, checked in by tnavarro, 14 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging
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
67      do ig=1,ngridmx
68
69         !write(*,*) "all qsurf to zero. dirty."
70         do iq=1,nqmx
71         qsurf(ig,iq)=0.  !! on jette les inputs GCM
72                          !! on regle juste watercaptag
73                          !! il faudrait garder les inputs GCM
74                          !! si elles sont consequentes
75         enddo
76         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
77     .        ( albedodat(ig) .ge. 0.26   ) )  then
78                 write(*,*)"outlier ",ig
79                 watercaptag(ig)  = .true.
80                 dryness(ig)      = 1.
81                 albedodat(ig)    = albedo_h2o_ice  !! pour output
82         else
83                 watercaptag(ig)  = .false.
84                 dryness(ig)      = 1.
85         endif
86         
87      endif
88#else
89       ! print*,'ngridmx',ngridmx
90
91      if (improvedicelocation) then
92     
93        if (ngridmx .eq. 738) then ! hopefully 32x24
94          print*,'water ice caps distribution for 32x24 resolution:'
95          watercaptag(1:9)    = .true. ! central cap - core
96          watercaptag(26:33)  = .true. ! central cap
97           
98        else if (ngridmx .eq. 3010) then ! hopefully 64x48
99          print*,'water ice caps distribution for 64x48 resolution:'
100          watercaptag(1:65)   = .true. ! central cap - core
101!          watercaptag(85)     = .true. ! central cap
102          watercaptag(83:85)  = .true. ! central cap - BIGGER
103          watercaptag(93:114) = .true. ! central cap
104!          watercaptag(254)    = .true. ! Korolev crater
105!          watercaptag(255)    = .true. ! Korolev crater --DECALE
106          watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG
107!--------------------------------------------------------------------         
108!          watercaptag(136)    = .true. ! outlier, lat = 78.75
109!          watercaptag(138)    = .true. ! outlier, lat = 78.75
110!          watercaptag(140)    = .true. ! outlier, lat = 78.75
111!          watercaptag(142)    = .true. ! outlier, lat = 78.75
112!          watercaptag(183)    = .true. ! outlier, lat = 78.75
113!          watercaptag(185)    = .true. ! outlier, lat = 78.75
114!          watercaptag(187)    = .true. ! outlier, lat = 78.75
115!          watercaptag(189)    = .true. ! outlier, lat = 78.75
116!          watercaptag(191)    = .true. ! outlier, lat = 78.75
117!          watercaptag(193)    = .true. ! outlier, lat = 78.75
118!--------------------------------------------------------------------         
119          watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75
120          watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75
121           
122        else if (ngridmx .ne. 1) then
123      print*,'No improved ice location for this resolution!'
124      print*,'Please set improvedicelocation to false in surfini.'
125      print*,'And check lat and lon conditions for ice caps in code.'
126      call exit(1)
127         
128        endif
129       
130        ! print caps locations
131        print*,'latitude | longitude | ig'
132        do ig=1,ngridmx
133          dryness (ig) = icedryness
134          !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs
135          !!! OU remettre comme c'était avant Aymeric et pis c'est tout
136          if (watercaptag(ig)) then
137             print*,'ice water cap', lati(ig)*180./pi,
138     .              long(ig)*180./pi, ig
139          endif
140        enddo
141
142!!------------------------ test -----------------------------
143!!-----------------------------------------------------------   
144!!      else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA ....
145!!       
146!!         do ig=1,ngridmx
147!!          dryness (ig) = 1
148!!          if (albedodat(ig) .ge. 0.35) then
149!!             watercaptag(ig)    = .true.
150!!             print*,' albedo criteria ice water cap', lati(ig)*180./pi,
151!!     .              long(ig)*180./pi, ig
152!!          endif
153!!        enddo
154!!-----------------------------------------------------------   
155!!-----------------------------------------------------------   
156     
157      else
158
159         do ig=1,ngridmx
160
161         dryness (ig) = icedryness
162         
163!!c Towards olympia planitia water caps ('relatively' low latitude ones)
164!!c---------------- proposition par AS --------------------
165!!c--------------------------------------------------------
166!!c          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
167!!c     .         ( lati(ig)*180./pi      .le.  77 ) .and.
168!!c     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
169!!c     .             ( long(ig)*180./pi .le. 120. ) )
170!!c     .           .or.
171!!c     .           ( ( long(ig)*180./pi .ge. -130. ) .and.
172!!c     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
173!!c---------------- proposition par TN --------------------
174!!c---------------- HIGHLY EXPERIMENTAL -------------------
175!!c--------------------------------------------------------     
176!!       if ( ( ( lati(ig)*180./pi .ge. 73.  ) .and. ! cap #1
177!!     .           ( lati(ig)*180./pi .le. 75.1 ) .and.
178!!     .           ( long(ig)*180./pi .ge. 95.  ) .and.
179!!     .           ( long(ig)*180./pi .le. 110. ) )
180!!     .         .or.
181!!     .         ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
182!!     .           ( lati(ig)*180./pi .le. 80.  ) .and.
183!!     .           ( long(ig)*180./pi .ge. 110. ) .and.
184!!     .           ( long(ig)*180./pi .le. 140. ) )
185!!     .         .or.
186!!     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
187!!     .           ( lati(ig)*180./pi .le. 78.  ) .and.
188!!     .           ( long(ig)*180./pi .ge. 155. ) .and.
189!!     .           ( long(ig)*180./pi .le. 180. ) )
190!!     .                .or.
191!!     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
192!!     .           ( lati(ig)*180./pi .le. 72.  ) .and.
193!!     .           ( long(ig)*180./pi .ge. 163. ) .and.
194!!     .           ( long(ig)*180./pi .le. 168. ) )
195!!     .         .or.
196!!     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
197!!     .           ( lati(ig)*180./pi .le. 78.  ) .and.
198!!     .           ( long(ig)*180./pi .ge. -160.) .and.
199!!     .           ( long(ig)*180./pi .le. -120.) ) )
200!!     .         then
201       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
202     .           ( lati(ig)*180./pi .le. 80.  ) .and.
203     .           ( long(ig)*180./pi .ge. 110. ) .and.
204     .           ( long(ig)*180./pi .le. 181. ) )
205     .         .or.
206
207     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
208     .           ( lati(ig)*180./pi .le. 76.  ) .and.
209     .           ( long(ig)*180./pi .ge. 150. ) .and.
210     .           ( long(ig)*180./pi .le. 168. ) )
211     .         .or.
212     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
213     .           ( lati(ig)*180./pi .le. 80.  ) .and.
214     .           ( long(ig)*180./pi .ge. -150.) .and.
215     .           ( long(ig)*180./pi .le. -110.) ) )
216     .         then
217     
218             if (temptag) then
219             
220               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
221                if (ngridmx.ne.1) watercaptag(ig)=.true.
222!                  print*,'ice water cap', lati(ig)*180./pi,
223!     .              long(ig)*180./pi, ig
224                  !dryness(ig) = 1.
225                  alternate = 1
226                else
227                  alternate = 0
228                endif !end if alternate = 0
229               
230             else
231             
232!               if (ngridmx.ne.1) watercaptag(ig)=.true.
233!                  write(*,*) "outliers ", lati(ig)*180./pi,
234!     .              long(ig)*180./pi
235     
236             endif ! end if temptag
237             
238           endif
239
240
241c Opposite olympia planitia water cap
242c---------------- proposition par AS --------------------
243c--------------------------------------------------------
244c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
245c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
246c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
247c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
248c---------------- proposition par TN --------------------
249c--------------------------------------------------------
250           if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
251     .            ( lati(ig)*180./pi     .le.  84 ) )
252     .          .and.
253     .          ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
254     .            ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
255!     .        ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
256!     .            ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
257!     .          ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
258!     .            ( long(ig)*180./pi .le. -70. ) ) ) ) then     !!! 64x48
259              if (ngridmx.ne.1) then
260                watercaptag(ig)=.true.
261                 print*,'ice water cap', lati(ig)*180./pi,
262     .            long(ig)*180./pi, ig
263              endif
264              !dryness(ig) = 1.
265           endif
266
267c Central cap
268c---------------- anciens reglages FF -------------------
269c--------------------------------------------------------
270
271           if (lati(ig)*180./pi.gt.84) then
272             print*,'ice water cap', lati(ig)*180./pi,
273     .         long(ig)*180./pi, ig
274             if (ngridmx.ne.1) watercaptag(ig)=.true.
275             !dryness(ig) = 1.
276c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
277c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
278c               if (ngridmx.ne.1) watercaptag(ig)=.true.
279c               dryness(ig) = 1.
280c             endif
281c             if (lati(ig)*180./pi.ge.85) then
282c               if (ngridmx.ne.1) watercaptag(ig)=.true.
283c               dryness(ig) = 1.
284c             endif
285           endif  ! (lati>80 deg)
286           
287         end do ! (ngridmx)
288         
289       endif ! of if (improvedicelocation)
290#endif     
291
292       ENDIF ! (caps & water)
293
294c ===============================================================
295c      INITIAL ALBEDO
296c ===============================================================
297
298         write(*,*)"surfini: water frost thickness",
299     s     frost_albedo_threshold
300         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
301         write(*,*)"surfini: water ice TI:", inert_h2o_ice
302
303c        To start with : Initial albedo = observed dataset
304c        -------------------------------------------------
305         DO ig=1,ngrid
306              psolaralb(ig,1)=albedodat(ig)
307              psolaralb(ig,2)=albedodat(ig)
308         END DO
309         PRINT*,'minimum albedo sans water caps',
310     s     albedodat(ISMIN(ngrid,albedodat,1))
311         PRINT*,'maximum albedo sans water caps',
312     s     albedodat(ISMAX(ngrid,albedodat,1))
313
314c        initial albedo if permanent H2O ice is present
315c        ------------------------------------------------
316         IF ((water) .and. (caps)) THEN
317           DO ig=1,ngrid
318            IF (watercaptag(ig)) THEN
319              psolaralb(ig,1) = albedo_h2o_ice
320              psolaralb(ig,2) = albedo_h2o_ice
321            ENDIF
322           END DO
323           PRINT*,'minimum albedo avec water caps',
324     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
325           PRINT*,'maximum albedo avec water caps',
326     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
327
328         ENDIF
329
330c      changing initial albedo if CO2 ice is present
331c      -------------------------------------------
332
333       DO ig=1,ngrid
334         IF (piceco2(ig) .GT. 0.) THEN
335             IF(ig.GT.ngrid/2+1) THEN
336                icap=2
337             ELSE
338                icap=1
339             ENDIF
340             psolaralb(ig,1) = albedice(icap)
341             psolaralb(ig,2) = albedice(icap)   
342         END IF
343       END DO
344
345c      changing initial albedo if water ice frost is present
346c      -------------------------------------------
347       IF (water) THEN
348          do iq=1,nqmx
349c          if there is frost and surface albedo is set to albedo_h2o_ice
350           if(noms(iq).eq."h2o_ice") then
351             do ig=1,ngrid
352              if ((piceco2(ig) .eq. 0.).and.
353     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
354                     psolaralb(ig,1) = albedo_h2o_ice
355                     psolaralb(ig,2) = albedo_h2o_ice
356c                     PRINT*,'surfini.F frost',
357c     &                  lati(ig)*180./pi, long(ig)*180./pi
358               endif
359              enddo
360           endif
361          end do
362          PRINT*,'minimum albedo avec givre et co2',
363     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
364          PRINT*,'maximum albedo avec givre et co2',
365     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
366       END IF
367         
368
369      RETURN
370      END
Note: See TracBrowser for help on using the repository browser.